Assessing the impact of large volcanic eruptions of the Last Millennium on Australian rainfall regimes

Explosive volcanism is an important natural climate forcing, impacting global surface temperatures and regional precipitation. Although previous studies have investigated aspects of the impact of tropical volcanism on various ocean-atmosphere systems and regional climate regimes, volcanic eruptions remain a poorly understood climate forcing and climatic responses are not well constrained. In this study, volcanic eruptions are explored in particular 15 reference to Australian precipitation, and both the Indian Ocean Dipole (IOD) and El Nino-Southern Oscillation (ENSO). Using nine realisations of the Last Millennium (LM) with different time-evolving forcing combinations, from the NASA GISS ModelE2-R, the impact of the 6 largest tropical volcanic eruptions of this period are investigated. Overall, we find that volcanic aerosol forcing increased the likelihood of El Nino and positive IOD conditions for up to four years following an eruption, and resulted in positive precipitation anomalies over northwest (NW) and southeast 20 (SE) Australia. Larger atmospheric sulfate loading during larger volcanic eruptions coincided with more persistent positive IOD and El Nino conditions, enhanced positive precipitation anomalies over NW Australia, and dampened precipitation anomalies over SE Australia.


Introduction
Volcanic eruptions have significant impacts on weather and climate variability through the injection of volcanogenic material into the atmosphere.Sulfate aerosols, formed through the reaction of SO 2 and OH -in the volcanic cloud, decrease incoming shortwave radiation, and if injected into the stratosphere, can generate a global response (Driscoll et 30 al., 2012;LeGrande et al., 2016).Previous studies have identified relationships between volcanism and surface and tropospheric cooling (Driscoll et al., 2012), local stratospheric warming (Wielicki et al., 2002), strengthening of the Arctic Oscillation and Atlantic meridional overturning circulation (Oman et al., 2005;Stenchikov et al., 2006, 2009& Shindell et al., 2004), and negative global precipitation anomalies (Gillet et al., 2004).The present study focuses on the under-studied relationship between large, globally significant tropical eruptions in the Last Millennium (850-1850CE) 35 and the El-Nino Southern Oscillation (ENSO), Indian Ocean Dipole (IOD) and Australian precipitation.
Previous studies have focused on the response of ENSO to volcanic eruptions.Adams et al. (2003) demonstrated a statistical relationship between explosive tropical volcanism and the ENSO; a tropical mode that has been linked to variations in Australian precipitation anomalies (Meyers et al., 2007;Pepler et al., 2014).They showed that large 40 tropical eruptions can increase the likelihood and amplitude of an El Nino event in following years, followed by a weaker La Nina state.Further work by Mann et al. (2005), Emile-Geay et al. (2008), McGregor et al. (2010), Wahl et al. (2014) and Predybaylo et al. (2017) supported this result.Pausata et al. (2015) identified that a radiative forcing threshold value of more than 15 W m -2 is required to affect the ENSO, and that high latitude Northern Hemisphere eruptions, in addition to tropical eruptions, are capable of doing so, as long as the forcing is asymmetric with regards to 45 the equator.
The relationship between volcanic forcing and ENSO has been attributed to two contrasting, though not unrelated, mechanisms.The dynamical thermostat mechanism (Clement et al., 1996), whereby a uniform reduction of the surface heat flux due to volcanism causes warming of the eastern equatorial Pacific, was identified as the driver of ENSO's 50 response to volcanism by Mann et al. (2005) and Emile-Geay et al. (2008).Conversely, a shift in the Intertropical Convergence Zone (ITCZ) induced by strong radiative forcing, was accredited in more recent studies (Pausata et al., 2015;Stevenson et al., 2016).Preconditioning does impact the severity of the ENSO response.Predybaylo et al. (2017) found that years with an initial central Pacific El Nino ENSO phase show the largest statistical impact from Pinatubosized eruptions and that summer eruptions coincided with a more pronounced El Nino response.

55
Despite the understanding that volcanism can trigger or amplify El Nino events in the following years, the exact relationship between ENSO and volcanic forcing is still debated.McGregor andTimmermann (2011) andZanchettin et al. (2012) reported an enhanced probability of La Nina events occurring in the immediate years after a volcanic eruption, rather than El Nino, while several other studies (Self et al., 1997;Robock, 2000;Ding et al., 2014) found no 60 relationship between ENSO and volcanic forcing.Robock (2000) argued that both El Chichon and Pinatubo reached their peak forcing after the initiation of El Nino events, indicating a coincidental relationship, while other studies (Driscoll et al., 2012;Lewis & Karoly, 2014;Lewis & LeGrande, 2015;Predybaylo et al., 2017) have pointed out challenges in determining long-term characteristics of ENSO due to short instrumental records, and its relationship to volcanic forcing due to variable representations of both ENSO and volcanic aerosols in GCMs (Global Climate 65 Models).
Comparatively little research has gone into the effects of volcanic forcing on the Indian Ocean Dipole (IOD), despite its known climatic impacts on Indian Ocean basin countries, such as Australia, South Africa, India and Indonesia (Cheung & Abram, 2016).The IOD is the zonal SST gradient between the tropical western Indian Ocean (WIO) and the tropical 70 south eastern Indian Ocean (EIO) (Roxy et al., 2011), defined by the Dipole Mode Index (DMI).Cheung & Abram (2016) found that the DMI shows a statistically significant correlation to volcanic forcing, with a negative IOD (nIOD) occurring immediately after an eruption and a positive IOD (pIOD) one year later.Maher et al. (2015) found a similar relationship, with coinciding El Nino and pIOD events occurring 6-12 months after the peak of volcanic forcing.The response of the IOD to volcanic forcing has been hypothesised to result from either the IOD's relationship with ENSO 75 (Cheung & Abram, 2016), or the volcanically-induced reduction of the Asian Monsoon (Anchukaitis et al., 2010;Zambri et al., 2017).
Australian precipitation is affected by both the IOD and ENSO.Positive IOD (pIOD) and El Nino events typically cause averaged precipitation deficits, while negative IOD (nIOD) events and La Nina cause positive precipitation 80 anomalies (Meyers et al., 2007;Pepler et al., 2014).While ENSO is often held responsible for triggering Australian droughts, the IOD has been shown to have an equal, if not larger, impact on heavily populated areas of Australia, with all significant southeastern Australia droughts in the 20th C showing a larger response to pIOD events than El Ninos (Ummenhofer et al., 2009).Volcanic aerosols have been found to cause global precipitation deficits for up to 5 years post-eruption (Robock & Lui, 1994;Iles et al., 2013;Gillett et al., 2004;Gu & Adler, 2011;Soden et al., 2002;Joseph & Zeng, 2011;Schneider et al., 2009;Timmreck et al., 2012;Iles et al., 2015).However, these deficits have been shown to vary seasonally (Joseph & Zeng, 2011), and cause positive precipitation anomalies over the NW and SE of Australia in the Southern Hemisphere (SH) winter and early spring (June-September), despite significant precipitation deficits in the summer (Joseph & Zeng, 90 2011;Schneider et al., 2009).Schneider et al. (2009) found a precipitation increase in the SH winter in Australia, Africa and Southern Asia.This current study explores these unclear relationships between volcanic eruptions and Australian rainfall.

Simulations
To understand the response of the IOD, ENSO and Australian precipitation towards volcanic forcing in the last millennium, we analysed 9 ensembles from the NASA GISS ModelE2-R (hereafter simply GISS) (Schmidt et al., 100 2014).The GISS ensemble was run for the pre-industrial part of the last millennium (LM), from 850-1850 CE, which is defined by the Coupled Model Intercomparison Project Phase 5 (CMIP5) (Taylor et al., 2012).GISS is run at 2 degrees x 2.5 degrees horizontal resolution, with 40 vertical levels up to 0.1 hPa.
Five ensembles were forced with volcanic forcing, while four were not.Of the five run with volcanic forcing, four were 105 forced with Crowley and Unterman (2013)'s aerosol optical depth data (CR), and one with double the Ice-core Volcanic Index 2 by Gao et al. (2008) (2xG) (see Table 1 for experiment summary).The GISS model is forced with prescribed Aerosol Optical Depth (AOD) from 15-35 km.The LM simulations also include transient solar and land use histories that differ between ensembles.However, as this analysis focuses primarily on the immediate post-volcanic response, the impact of these smaller amplitude and slowly varying forcings is likely to be insignificant (Colose et al., 2016).110

Methods
First, the six largest tropical eruptions between 850-1850CE were identified by the magnitude of their total global stratospheric sulfate aerosol injection (Tg) from the IVI2 Version 2 dataset, revised in 2012 (Gao et al., 2008), and the 115 years surrounding eruption extracted for analysis (see Table 2 and Figure 1).Eruptions were deemed as tropical if volcanic aerosols were present in significant amounts in both hemispheres.
The Kuwae eruption is included within the analysed eruptions, and is dated to 1452CE.While this year contains the bihemispheric deposition from the Kuwae eruption in both volcanic datasets used in this study, it is important to note that 120 Sigl et al. (2013) recently constructed an ice-core record of volcanism that dates the Kuwae eruption to 1459/1459CE, with another, smaller eruption occurring at 1452.For the purposes of this paper, however, Kuwae will be considered as the 1452 deposition event.
We explored anomalous conditions in ENSO, the IOD and Australian rainfall.For ENSO, the period December-125 February (DJF) was examined using the NINO3.4index, defined by the averaged sea surface temperature (SST) anomalies between 5N-5S and 170-120W.When analysing the IOD, the July-November (JASON) period was examined due to the IOD's tendency to develop and mature over these months (Weller et al., 2014).The IOD was measured using  Australian precipitation was processed to find the anomalies of each season and year relative to the long-term mean.
Rainfall was also examined over the period JASON.Analyses were conducted on south-eastern Australia (132.5-155E; 27.5-45S) and north-western Australia (110-132.5E;10-27.5S).The south-east and north-west were chosen for analysis as the effect of the IOD on Australian precipitation is largest in, and potentially limited to, these areas (Ashok et al., 135 2004).
The response of these large-scale modes of variability and rainfall are investigated using an epoch approach.For each major identified eruption, a response was defined by subtracting a reference period (the mean of 5 years pre-eruption) from the eruption year and the six years following eruption individually.A reference period of five years was chosen as 140 it minimised the effect of trends or low-frequency climate variability (Iles et al., 2013).Mode specific graphs (IOD, ENSO, Australian Precipitation) focused on the nine years surrounding eruption (years -2 to 6, with year 0 being the year of eruption).The mean of all six eruptions in each ensemble were calculated for individual years, and then the mean of all ensembles included in each forcing category were compared (CR, 2xG or None).145

Results
The global SST response for the CR forcing group shows predominantly surface cooling anomalies (Fig. 2).More specifically, cooling occurs in the Northern Hemisphere from years 0-3, while the south Atlantic, Indian and Pacific oceans show mostly minor warming.An El Nino-like temperature gradient is apparent in years 0,1 and 4, relative to the 150 eruption, despite overall cooling trends.In the Southern Hemisphere, cooling is most pronounced over land masses, particularly Australia and the southern tip of Africa.
The DMI response showed a significant pIOD condition one year after a major eruption in all volcanically forced ensembles that persists until year 5, where an abrupt negative IOD phase occurs (Fig. 3).This response can also be seen 155 in Fig. 2, where the EIO region shows larger and more widespread cooling anomalies than the WIO region in years 1, 2 and 4.This response contrasts to the non-volcanically forced ensembles, which show neither a prolonged pIOD nor nIOD condition.The response of the DMI to the largest and smallest eruptions were also extracted.Fig. 4 show the mean DMI response to the 1258 Samalas eruption (257.91 Tg) and the 1600 Huaynaptina eruption (56.59 Tg).Our results show that while both eruptions caused a significant simulated pIOD at year 1, the larger 1258 Samalas eruption 160 persisted with a significant pIOD condition in years 2 and 4, while the 1600 Huaynaptina eruption did not.
The mean NINO3.4 multi-volcano response to ensemble forcing showed a statistically significant El-Nino like response for all 6 years following eruption, with a peak at year 3 in both the CR and 2xG ensembles (Fig. 5).The nonvolcanically forced ensemble group showed neither a significant El Nino nor La Nina tendency, with the NINO3.4165 index remaining within 0.4/-0.4.The index also showed an increase in the intensity and endurance, of post-volcanic El Ninos between the Samalas and Huaynaptina eruptions (Fig. 6).The Samalas eruption was followed by an El Nino that endured for 3 years, from years 1-3, peaking at a NINO3.4anomaly of 0.68 in year 3, while Huaynaptina peaked at 0.53 in year 2 from an El Nino that endured for 2 years.Ensembles with volcanic forcing showed an increase in precipitation over southeast (SE) (Fig. 10) and northwestern (NW) (Fig. 8) Australia between July to November (JASON).Both areas showed predominantly positive anomalies in 180 years 0-5 post-eruption, with the largest response seen between years 0-2.NW Australia (Fig. 8) showed larger positive precipitation anomalies between years 0-2 than SE Australia (Fig. 10) in the CR ensemble mean, and in years 0 and 2 in the 2xG ensemble mean.
Comparison of the precipitation anomalies following the Samalas and Huaynaptina eruptions in NW Australia (Fig. 9) 185 showed that the smaller eruption had a delayed and smaller positive precipitation peak, with Samalas peaking in year 0 with an anomaly of 0.23 and Huaynaptina in year 2 at 0.14.While the Huaynaptina eruption also showed a delayed peak in precipitation in SE Australia (Fig 11), the persistence of positive precipitation anomalies exceeded those of the Samalas eruption.Huaynaptina recorded values >0.17 in years 1-2 and a value of 0.12 in year 4, all of which were larger anomalies than the peak of the Samalas eruption at 0.11 in year 1 (Fig. 11).190

Discussion and conclusions
Our results suggest that the large-scale IOD and ENSO systems, and Australian rainfall regimes, were all impacted by large tropical eruptions of the Last Millennium.

195
The DMI response simulated in the GISS ensemble following large eruptions is complimentary to previous research conducted by Cheung & Abram (2015) and Maher et al. (2015).The pIOD peak in year 1 (Fig. 3) is consistent with both studies, in which statistically significant pIOD conditions occurred between 6 months to 2 years after an eruption.Cheung & Abram (2015) also found a statistically significant negative condition immediately after eruption at year 0, 200 however this was absent from both Maher et al.'s (2015) results and the CR forcing category in this study.The 2xG category does show a nIOD condition at year 0, but is not believed to be a response to volcanic forcing as a similar nIOD condition can be seen at year -1.The abrupt shift to a negative condition at year 5 was not found in either Cheung & Abram (2015) or Maher et al. (2015)'s results.Both studies found a gradual decrease in DMI from year 1 to years 3-4.

205
The smooth transition to a lower DMI following eruptions found by Cheung & Abram (2015) and Maher et al. (2015) contrasts with the abrupt change from a pIOD of approximately 0.13 at year 4, to an nIOD of -0.069 in the CR ensembles and -0.083 in the 2xG ensemble at year 5 (Fig. 3).This inconsistency between studies could be due to the selection of eruptions analysed by each paper.Cheung & Abram (2015) included all eruptions from 850-2005CE 210 recorded on the IVI2 in their analysis.While this encompasses all eruptions analysed here, it also included many smaller eruptions that would likely have dampened the climatic response, a response that has been analysed in previous papers (Zambri & Robock, 2016).Maher et al. (2015) looked at the five largest eruptions from 1880 to present, of which the largest was Pinatubo (1991), measured at 30.10 Tg globally on the IVI2 (Gao et al., 2008).In comparison, our research deals with eruptions of much larger atmospheric loading, ranging from 56.59 to 257.91 Tg.Therefore, the 215 smallest eruption dealt with in this paper is almost double that of the largest in Maher et al.'s (2015) analysis.
The persistence of a high pIOD through to year 5 seen here may result from the larger mean atmospheric sulfate loading imposed.This theory is supported by the comparison between the Samalas and Huaynaptina (Fig. 4) eruptions.Our results showed that while both eruptions caused a significant pIOD at year 1, the larger 1258 Samalas eruption alone 220 persisted with a significant pIOD condition in following years.Further support can be gathered from the comparison between the 2xG and CR ensemble means in fig. 3. Years 0-3 show more extreme values in the 2xG ensemble mean, while years 4-6 show similar values for both forcing categories.Maher et al. (2015) found a similar response, with the two largest eruptions analysed in the paper showing the largest and longest enduring pIOD anomalies.This suggests that larger mean atmospheric sulfate loading can cause not only more persistent, but also more extreme pIOD  2010) and Maher et al. (2015).Fig. 5 shows a very prominent and persistent El-Nino response in all 6 years following eruption, however it lacks the weaker La Nina-like state that was observed 3-6 230 years after eruption in these previous papers.Spatial maps of our SST (Fig. 2), while dominated by the overall volcanic cooling, show an El Nino-like pattern in the eastern Pacific that is most visible in year 4, but is also distinctive in years 0, 1 and 3. We can therefore conclude that an El Nino-like anomaly was generated in the multimodel mean response in years 0-6 following eruption, despite the overall cooling trend caused by the volcanic aerosols.Comparison of the Samalas and Huaynaptina (Fig. 6) eruptions also suggest that, similar to the DMI, the intensity and endurance of the 235 ENSO response to volcanic forcing increases with increasing mean atmospheric sulfate loading.This once again supports the findings of Maher et al. (2015) that identified a similar pattern.
The positive response of Australian precipitation to volcanic forcing as seen here (Fig. 7, 8 &10) is in agreement with several papers that identified positive precipitation responses over Australia to large volcanic eruptions (Schneider et 240 al., 2009;Joseph & Zeng, 2011).Our results suggest that the direct precipitation effect of volcanic aerosols override the impact of the IOD on Australian precipitation in the years following large tropical volcanic eruptions.NW Australia (Fig. 8) showed larger positive precipitation anomalies between years 0-2 than SE Australia (Fig. 10) in the CR ensemble mean, and in years 0 and 2 in the 2xG ensemble mean.This could be due to the positive precipitation anomalies that can be generated by combined El Nino and pIOD events in the NW Australian region (Meyers et al., 245 2007& Pepler et al., 2014), enhancing the precipitation surplus caused by volcanic aerosols.
The varying response of NW Australia to the Samalas and Huaynaptina eruptions (Fig. 9) also supports the enhancement of the volcanically induced precipitation surplus by combined El Nino and pIOD events.The Samalas eruption was followed by strong and enduring El Nino and pIOD conditions for up to 4 years post volcanism, and 250 showed larger positive precipitation anomalies from years 0-3 than the Huaynaptina eruption, that was accompanied by smaller, shorter-lived El Nino and pIOD conditions.The precipitation surplus to the Samalas eruption in NW and SE Australia also peaked earlier than Huaynaptina, which could be a response to the larger atmospheric sulfate loading.425 Zanchettin, D., Timmreck, C., Graf, H. F., Rubino, A., Lorenz, S., Lohmann, K., Krüger, K., and Jungclaus, J. H.: Bi-decadal variability excited in the coupled ocean-atmosphere system by strong tropical volcanic eruptions, Clim. Dyn., 39, 419-444, 2012 Clim.Past Discuss., https://doi.org/10.5194/cp-2017-109Manuscript under review for journal Clim.Past Discussion started: 25 September 2017 c Author(s) 2017.CC BY 4.0 License.

Fig. 7
Fig.7shows the mean precipitation response of all volcanically forced ensembles.Substantial precipitation deficits can be seen in the tropics in years 0-2, with a deficit band of decreased precipitation also occurring at approximately 40°S and between 0-40°N in the western Pacific and North Atlantic Oceans.The Southern Hemisphere subtropics appear to have a slight increase in precipitation in years 0-2, most prominently, over Australia and southernmost Africa, while the southern polar region (60-90°S) shows only variable minor precipitation anomalies occurring in all 6 years post 175 Clim.Past Discuss., https://doi.org/10.5194/cp-2017-109Manuscript under review for journal Clim.Past Discussion started: 25 September 2017 c Author(s) 2017.CC BY 4.0 License.

Fig 7 :
Fig 7: Global precipitation anomalies (mm/day) in response to the Crowley and Unterman (CR) forcings, showing multimodel mean responses averaged across all analysed eruptions for years 0 to +5 after eruption.

Fig 9 :Fig 10 :
Fig 9: Mean NW Australian precipitation (mm/day) response across all volcanic ensembles to the largest (1258 Samalas) and smallest (1600 Huaynaputina) eruptions analysed.The bold lines represent the mean of all volcanic ensembles to each eruption, and the fainter lines represent the 90 th and 10 th percentile.

Fig 11 :
Fig 11: Mean SE Australian precipitation (mm/day) response across all volcanic ensembles to the largest (1258 Samalas) and smallest (1600 Huaynaputina) eruptions analysed.The bold lines represent the mean of all volcanic ensembles to each eruption, and the fainter lines represent the 90 th and 10 th percentile.495