The sulfur- and halogen-rich super eruption Los Chocoyos and its impacts on climate and environment

. The super-eruption of Los Chocoyos, newly dated to 80.6 kyrs ago, in Guatemala was one of the largest volcanic events of the past 100,000 years. Recent petrologic data show that the eruption released very large amounts of climate-relevant sulfur and ozone destroying chlorine and bromine gases. Using the recently released Earth System Model CESM2(WACCM6) we simulate the impacts of the sulfur- and halogen-rich Los Chocoyos (~15° N) eruption on the preindustrial Earth System for the eruption month January. Our model results show that enhanced modeled sulfate burden and aerosol optical depth (AOD) persists for five years, while the volcanic halogens stay elevated for nearly 15 years. As a consequence the eruption leads to a collapse of the ozone layer with global mean column ozone values dropping to 50 DU (80 % decrease) leading to a 550 % increase in surface UV over the first five years with potential impacts on the biosphere. The volcanic eruption shows an asymmetric hemispheric response with enhanced aerosol, ozone, UV, and climate signals over the Northern Hemisphere (NH). Surface climate is impacted globally due to peak AOD of > 6 leading to a maximum surface cooling of >6 K, precipitation and terrestrial net primary production (NPP) decreases of >25 %, and sea ice area increases of 40 % in the first three years. Locally, a wetting (>100 %) and strong increase of NPP (>700 %) over Northern Africa is simulated in the first five years related to

model uncertainties for the climate response to super eruptions are very large observational evidence from paleo archives and a coordinated model intercomparison would help to improve our understanding of the climate and environment response.

Introduction
The Los Chocoyos (LCY) super-eruption (Magnitude Mv=8 (Kutterolf et al., 2016), dated to 80.8±6.7 kyrs before present (Cisneros et al. to be submitted), was already 30 years ago assumed to be one of the largest volcanic eruptions of the past 100,000 years (Drexler et al., 1980). Originating from present-day Guatemala, the eruption formed the current stage of the large Atitlán caldera. Los Chocoyos released more than ~1100 km 3 of tephra and the eruption is used as a widespread key stratigraphic marker during that time (Cisneros et al., to be submitted; Kutterolf et al., 2016). The ash layers can be found in marine deposits from offshore Ecuador to Florida over an area of more than 10 7 km2 (Kutterolf et al., 2016). Hardly anything is known about the climate impacts of this eruption from proxy records, but LCY emitted large amounts of climate and environmentally relevant gases including sulfur, chlorine and bromine compounds Kutterolf et al., 2015Kutterolf et al., , 2016Metzner et al., 2014).
The sulfur gases emitted by volcanoes have a strong direct climate impact through the formation of sulfuric acid aerosols which blocks incoming sunlight and cools the surface (Robock, 2000). Halogen compounds, like chlorine and bromine, contributes to catalytic ozone depletion in the stratosphere (Brasseur and Solomon, 2005;Solomon, 1999). There is welldocumented petrological evidence that large to extremely large explosive volcanic eruptions emitted environmentally significant amounts of chlorine and bromine (Cadoux et al., 2015(Cadoux et al., , 2018Krüger et al., 2015;Kutterolf et al., 2013Kutterolf et al., , 2015Vidal et al., 2016). Furthermore, recent atmospheric observations revealed that even relatively small volcanic eruptions can inject significant amounts of halogen compounds in the stratosphere (for review and overview discussions see (von Glasow et al., 2009;Krüger et al., 2015;WMO, 2018). Next, a sulfur and halogen rich eruption is expected to cool the Earths' surface and potentially damage the stratospheric ozone layer, with further impacts on the surface environment through the change in the atmosphere's transparency to harmful ultraviolet (UV), particularly shortwave UV-B, radiation (i.e. Brenna et al., 2019).
Even though there is little literature about LCY, the climate impact of super-eruptions has been discussed in the scientific literature since at least the early 1990s. Early studies argued that the eruption of Toba (73 kyrs) could have initiated a glacial period Self, 1992, 1993;Zielinski et al., 1996). In addition, there is evidence that human populations went through a genetic bottleneck (i.e., most of the population died) at approximately the same time as the eruption of Toba (Ambrose, 1998;Haslam and Petraglia, 2010;Williams et al., 2009), which has since been ruled out as unlikely (Timmreck et al., 2010. A thorough investigation of the climate impacts of very large volcanic eruptions (Mv: 7-8) requires the use of comprehensive coupled climate models or, ideally, Earth System Models (ESMs). There are several published simulation studies of such eruptions, mostly focusing on the Toba eruption, and the climate impact of it's sulfur injection. Bekki et al. (1996) used a two-dimensional chemistry transport model with internally generated atmospheric circulation to study the Toba eruption impact on the atmosphere. Their simulations indicate that Toba could have caused a long-lasting atmospheric response due to the interactions between chemistry and aerosol microphysics.
Later, Jones et al. (2005) used a coupled atmosphere-ocean general circulation model (AOGCM) to study the Toba eruption impact on climate. In this study they forced the model by linearly scaling the observed aerosol optical depth (AOD) from the 1991 eruption of Mt. Pinatubo by a factor of 100, resulting in peak cooling of ~11 K in the 2 nd post-eruption year, followed by an initial recovery taking a decade. Surface cooling larger than 1 K persisted for more than 20 years. The volcanic aerosol forcing only lasted 5 years, so the response needed to be maintained by feedbacks in other components of the Earth System through, i.e., sea-ice/ocean feedbacks sustaining the short atmospheric forcing to longer (decadal to centennial) time scales (Miller et al., 2012;Stenchikov et al., 2009;Zhong et al., 2011).
In a similar study, Robock et al. (2009) used three different AOGCMs to study the Toba eruption effect on climate and ozone. In their study, both the linear scaling AOD approach and directly injecting sulfur into a model with an interactive bulk aerosol module were utilized. The resulting magnitude and length of the cooling was similar to what was published in Jones et al. (2005) across the different model versions and forcing methods. The scenarios representing a 100 times Pinatubo forcing resulted in ~12 K peak cooling with multi-decadal recovery times. In a 300 times Pinatubo scenario simulation including atmospheric chemistry effects, Robock et al. (2009) found slightly stronger (~1 K) and much more long-lasting surface cooling (length of >10 K cooling extended by ~5 years) compared to the similar forcing scenario without chemistry, due to depletion of atmospheric hydroxyl (OH) limiting the speed of aerosol formation, leading to a longer-lasting forcing. In addition, they simulated an increase in global column ozone, attributed to the reduction of reactive hydrogen oxides in the atmosphere.
Another model study of Toba, presented in Timmreck et al. (2010, 2012, simulated a smaller climate impact with peak cooling of 3.5 K lasting up to 10 years from injected sulfur compared to the analogous simulations in Robock et al. (2009). A key difference between Timmreck et al. (2010) and previous AOGCM simulations of Toba is the inclusion of online aerosol microphysics in a modal aerosol scheme and the OH limitation mechanism, leading to larger aerosol sizes, lower peak AOD (~4) values and thus lower climate impacts in Timmreck et al. (2010). The inclusion of interactive OH chemistry in the formation of aerosol is important because the availability of OH controls the speed of SO2 oxidation into sulfate (Bekki, 1995).
Concentrating on the atmospheric processes, English et al. (2013) used a sectional aerosol microphysical model coupled to a chemistry climate model with prescribed sea surface temperatures (SSTs) to study the aerosols and atmospheric impacts of the Toba eruption. In their model setup, neglecting aerosol radiative heating, they simulate even lower peak AOD values (~2.6) due to the fact that their sectional aerosol module creates larger aerosols compared to Timmreck et al. (2010).
In the climate modeling literature on the Toba super eruption there is a progression from larger climate (and environmental) impacts to smaller as models complexity develop over time. One key seems to be that climate effects are self-limiting for larger eruptions due to an increase of aerosol growth which reduces peak AOD (English et al., 2013;Pinto et al., 1989;Timmreck et al., 2010). In addition, the role of atmospheric chemistry and OH limitation on sulfuric acid aerosols is continuously under discussion in the literature (Bekki, 1995;Mills et al., 2017;Niemeier et al., 2019;Robock et al., 2009;Timmreck et al., 2003).
Investigating the effects of the Toba eruption on the earth system such as hydrology and terrestrial net primary production reveals a substantial reduction of precipitation globally leading to reduction of tree cover, increase of grass cover and decreased NPP, but with large regional and inter-model variability (Robock et al., 2009;Timmreck et al., 2010Timmreck et al., , 2012. While tropical deciduous trees and broadleaf evergreen trees virtually disappear in the simulations of Robock et al. (2009), Timmreck et al. (2011 find much more muted impacts on tree cover, particularly in the Tropics. Timmreck et al. (2011) even simulate an increase of NPP in the tropical rain forest regions of South America and Africa.
We are not aware of studies of super-size eruption effects on sea ice/ocean and the El Niño Southern Oscillation (ENSO), whereas the effects of large to very large volcanic eruptions have been widely discussed in the literature. Recent studies proposed a sea ice/ocean mechanism which prolong the volcanic induced short, abrupt surface cooling and sea ice increase to longer time scales (decadal) with the ocean sustaining the signal by buffering and transporting the cooling poleward (Miller et al., 2012;Zhong et al., 2011). With regard to volcanoes and ENSO, there is an ongoing debate (Stevenson et al., 2017) that tropical volcanic eruptions can either lead to La Niña-like response in the same year (Anchukaitis et al., 2010;Li et al., 2013) or the following (five) years  or to El Niño-like response up to the two following years (e.g., (Adams et al., 2003;Handler, 1984;Khodri et al., 2017;Ohba et al., 2013;Stevenson et al., 2016)). Discussions include the significance and mechanism of the results as well as the eruption characteristics, latitude, season, and strength.
Simulations of halogen (and sulfur) rich eruptions show that these can have serious, long-lasting impacts on the ozone layer, with implications for the surface environment through the increase of surface ultraviolet radiation. In a study of the ~1 Myr period of extremely large effusive eruption of the Siberian Traps (~250 million years ago) (Black et al., 2014), a near-total elimination of the ozone layer was found, as well as massive increases in surface ultraviolet radiation (~10 years recovery time after cessation of emissions), Black et al. (2014) hypothesized that this could be a contributing mechanism to the hemisphere, with peak depletion of 20-90 % depending on the degassing budget. In another study using a 2D CTM (Klobas et al., 2017), (hypothetical) volcanic halogens were included with the sulfur injection of Mt. Pinatubo, showing ozone depletion of 20 % lasting a few years under different future emission scenarios. Using CESM1(WACCM) a comprehensive coupled chemistry climate model (CCM), with prescribed volcanic aerosols and SSTs, Brenna et al. (2019) simulated an average Central American Volcanic Arc (CAVA) eruption with magnitude 6.4. They found ozone depletion up to 20 % globally lasting up to 10 years, which were most pronounced over the Northern Hemisphere (NH) and dropped below present day ozone hole conditions over Antarctica and the tropics. Consequently, surface UV radiation increased by >80 % over the 2 years with potential impacts on human health, agriculture and marine life.
Despite the Siberian Trap volcanism studies (Beerling et al., 2007;Black et al., 2014), we are not aware of super eruption studies taking the combined effect of sulfur and halogen injections in a fully coupled aerosol chemistry climate/ earth system model into account.
In this study we use the recently released CESM2 coupled with WACCM6 as the atmospheric component, which allows us to investigate newly the coupling and the feedbacks between volcanic aerosols, chemistry, radiation, climate and the earth system after a sulfur and halogen rich super volcanic eruption. In this paper, we present the climate and environmental impacts of the sulfur-and halogen-rich super eruption Los Chocoyos (Guatemala, 80.6 kyrs ago; Cisneros et al. to be submitted). In a forthcoming paper, we will investigate the impacts on the atmospheric circulation.

Los Chocoyos erupted volatile estimates
Using the recently published total erupted mass estimate for the LCY eruption (Kutterolf et al., 2016) and the previously published petrologic estimates of volatile concentrations for sulfur, chlorine and bromine (Kutterolf et al., 2013Metzner et al., 2014) we calculate a new mass of erupted volatiles for LCY as a starting point for defining the stratospheric injections in our model simulations. The erupted volatile masses as calculated using these estimates are 523±94 Mt sulfur, 1200±156 Mt chlorine and 2±0.46 Mt bromine.

CESM2(WACCM)
In this study we use the Community Earth System Model Version 2 (CESM2) (Danabasoglu et al., (Liu et al., 2016), which has been adapted and extended for the stratosphere .
CESM2(WACCM6), as a coupled CCM, allows us to explore the coupling between radiation, temperature, circulation, chemistry and composition in the atmosphere. The horizontal resolution is 0.95° longitude by 1.25° latitude with 70 hybrid sigma pressure layers from the surface to 5.5e-6 hPa. The quasi-biennial oscillation (QBO) is internally generated

Model experiments
To model the impact of the LCY eruption on the atmosphere we use the petrological estimated erupted sulfur and halogen masses as input. We inject all of the erupted sulfur mass as SO2 and 10 % of the erupted halogen mass, as HCl and HBr, into the stratosphere, which we consider a conservative estimate given the range of available observational and plume model evidence (see the detailed discussion in Krüger et al. 2015 andBrenna et al. 2019). The volcanic volatiles are injected into the model grid boxes at 14.6°N and between 80° and 97.5°W, at 24 km altitude. The eruption date is set to 1. January, as the eruption season is not known, and needed to be spread over the first 6 days in the model to allow the massive volatile mass to accumulate over time and neighboring longitudinal grids.
We run an ensemble of six simulations for the combined sulfur and halogen forcings (LCY_full) starting from different ENSO (positive, negative, neutral) and QBO (easterly, westerly) states of the control run (Ctr). In addition, we perform two simulations (QBO easterly and westerly, ENSO neutral, branching from the same Ctr years as for LCY_full) with only the sulfur forcing (LCY_sulf) to explore the difference in response between the forcing scenarios. All ensemble members last 35 years. The control simulation is 70 years with constant 1850 forcings. In all simulations, background stratospheric concentrations of chlorine and bromine are 0.45 ppbv and 10.2 pptv at the 10 hPa level respectively, consistent with preindustrial estimates (WMO, 2014). In addition to these main experiments, we have performed two sensitivity simulations.
One with sulfur injection reduced by a factor of 100 (LCY_1%sulf) and one with full sulfur injection and 1 % halogen injection efficiency instead of 10 % (LCY_1%halog). The experiments are summarized in Table 1.

Surface UV calculations
We calculate the UV radiation at the surface using the Tropospheric Ultraviolet and Visible (TUV) radiation transport model (Madronich and Flocke, 1997)  hourly temporal resolution to get a representation of the variations in UV throughout the year. As input to the TUV model we give averages of column ozone and AOD over the first five years after the eruption for LCY_full and LCY_sulf and for the control run to generate UV fields for the eruption scenarios and for the control run climatology.

Oceanic Niño Index (ONI)
To select initial conditions for the set-up of the ensembles and to quantify the impact of the volcanic eruptions on the ENSO we calculate the Oceanic Niño Index (ONI) using the model output. An ENSO positive or negative phase is defined as ONI>0.5 or ONI<-0.5, respectively, for more than 5 consecutive months.

Atmospheric burdens of volcanic gases and aerosols
Using our modeling approach results in the atmospheric burdens of volcanic gases and aerosol distinguishing between burden anomalies and the same burdens, normalized by the respective maximum values as summarized in Figure 1.
The normalized total sulfur burden after the simulated eruptions has an e-folding time (reduction by 1/e) of a little more than two years (Figure 1b). There is first a plateau for ~1 year before decay starts. The following e-folding times (1/e 2 and 1/e 3 ) are shorter, a bit less than 1 year. After ~5 years most of the sulfur has been removed from the atmosphere. The total sulfur burden life times are remarkably similar for all four injection scenarios, even when the injection amounts differ by a factor of 100. In contrast, the halogens have a longer e-folding time of approximately 3 years, but the following e-folding times for 1/e 2 and 1/e 3 are ~1 year, similar to what we simulate for total sulfur.
The conversion of injected SO2 into sulfuric acid aerosol is significantly slowed down in the LCY_sulf, LCY_full and LCY_1%halog scenarios compared to the LCY_1%sulf, where we only injected 1 % of the sulfur (Figure 1 d). The e-folding time of the SO2 burden increases from ~3 months to ~6 months when the injected SO2 mass is increased by a factor of 100.
There is an additional increase in lifetime to ~1 year when halogens are injected in addition to SO2. This increase in SO2 lifetime is caused by the limited oxidizing capacity of the atmosphere (see also Bekki et al. 1995). The main compound responsible for oxidizing SO2 to sulfuric acid is OH. When the SO2 burden increases the availability of OH is limited and oxidation slows down. When halogens are injected in addition, reactions involving halogens also consume OH (Fig. 1f) further limiting the OH available for SO2 oxidation, and thus further increasing the SO2 lifetime (Fig. 1d).
The peak aerosol mass, when the sulfur burden is the same, is depending on the conversion time of SO2 into aerosol. Thus, the peak aerosol mass is lower in the LCY_full scenario compared to the LCY_sulf (Figure 1c (Figure 1e), as expected when the injection is smaller (English et al., 2013;Pinto et al., 1989;Timmreck et al., 2010). Even though the aerosols are smaller in the LCY_1%sulf simulation, the removal time scale for the aerosols is similar to the two other scenarios, This indicates that gravitational settling, which is not important for sub-micrometer aerosols (Seinfeld and Pandis, 2016), is playing a minor role as a removal mechanism for the aerosol mass, and removal processes will tend to happen on the transport time-scale of the stratosphere.
Since the maximum aerosol mass in the LCY_sulf is ~25 % larger compared to the maximum mass in LCY_full, while the aerosol sizes are approximately the same, we find that the peak AOD at 550 nm is much larger in the LCY_sulf scenario (>8) compared to the LCY_full scenario (~6) (Figure 2a). This translates into a larger energy imbalance at the top of the model in LCY_sulf (Figure 2b). The maximum radiative imbalance at the top of the model is approximately -50 W/m 2 in the LCY_sulf and -40 W/m 2 in the LCY_full scenario. In both cases, an initial strong negative net imbalance is followed by a small positive net imbalance after ~3.5 years, and throughout the climate recovery period.

Ozone and UV response
Global ozone collapses after the eruption in the LCY_full scenario, with whole column values decreasing by >80 % to a global mean value of 50 DU (a measure of the thickness of the ozone layer) during years 2-3 after the eruption (Figure 3a).
Ozone levels lower than present-day Antarctic ozone hole conditions (<220 DU) persist for 8 years over the whole globe ( Figure 3a, S1). Depletion shows a bi-modal distribution in the stratosphere, with maximum depletion in the upper (~4 hPa) and lower (~30 hPa) stratosphere ( Figure S2). In the lower stratosphere, where most of the ozone mass is located, >70 % of ozone is destroyed after 1 year, and this level of depletion persists for 7 years (not shown). Peak depletion in the lower stratosphere is >95 %. Significant global mean ozone column reduction lasts for ~12 years. In the Antarctic, ozone hole conditions continue re-occurring annually for 16 years ( Figure S1b). Compared to LCY_full, the ozone response in LCY_1%halog is smaller, but reveals a similar response. A substantial decrease to global mean column values of ~150 DU and a recovery after about 10 years; a larger and longer lasting ozone response as was simulated for an average CAVA eruption with magnitude 6 and 10 % halogen injection efficiency (Brenna et al. 2019).
In contrast, in the LCY_sulf scenario, the column ozone increases by more than 100 DU in the first year after the eruption ( Figure 3a). This is caused by the increase in heterogeneous chemistry taking place on the sulfate aerosols which reduces the concentrations of ozone-destroying NO x (Tie and Brasseur, 1995) and was modelled for large and super-size volcanic sulfur injections into a pre-industrial stratosphere with low chlorine background levels (Muthers et al., 2015;Robock et al., 2009).
The ozone increase decays in about 3 years and is only slightly elevated after that until post-eruption year 10 ( Figure 3a).
The increase in ozone is concentrated in the mid to high latitudes and mostly in the NH ( Figure S1c).
In Figure 4 we present global maps of total AOD (a,b), as well as anomalies and the climatologies for column ozone and in the extratropics compared to a band of low AOD in the Southern Hemisphere (SH) tropics, and the largest impacts in the NH (Figure 4 a, b). In LCY_full, AODs are smaller over Antarctica than at lower latitudes. This might be as transport to the Antarctic region is suppressed by the strengthening of the Westerlies winds surrounding the southern polar vortex (not shown here), which acts as a transport barrier for the volcanic aerosols (Toohey et al., 2013). In LCY_sulf the global mean AOD is larger (c.f. Figure 2a), which holds for the local distribution over the globe as well (Figure 4a+b). In the LCY_full scenario, the largest increases in surface UV are more than 1400 % in the Arctic and more than 1000 % in the Antarctic. Changes are generally smaller towards the equator, but no part of the planet experiences less than a 200 % increase. Global average increase over the five-year period is 545 % (Figure 3b). By contrast, in the LCY_sulf scenario the UV-B decreases by more than 80 % in the mid-to-high latitudes of the NH and by >60 % over most of the rest of the planet (Figure 4g,h). The UV response in our calculations are impacted by the ozone levels and the AOD (Figure 4a-

f), and in
LCY_full the AOD and ozone effects are opposing each other with the ozone effects being strongest, while in LCY_sulf they are both contributing to a decrease in surface UV.

Global surface temperature decreases for 30 years
Time series of global mean surface temperatures are shown in Figure 3b. For both scenarios, global mean surface temperature decreases more than 6 K and is back to climatological background after approximately 30 years. The difference between the LCY_full and LCY_sulf forcing scenarios is >1 K at peak cooling in year 2. If the aerosol response from the sulfur injection (which is the same in these two scenarios) was the same, we would expect the temperature response to be very similar, and perhaps LCY_full would be slightly cooler since ozone depletion in the stratosphere is a negative radiative forcing on the global climate system (Myhre et al., 2013). Instead, we interpret this difference in surface temperature response due to the large difference in peak AOD (Figure 2a).
In Figure 5 (a, b) we show pentadal average maps of surface temperature anomaly for both the LCY_full and LCY_sulf scenarios. Temperature anomaly patterns are relatively similar between the scenarios with surface cooling almost globally and largest anomalies in the NH and over the continents. The magnitudes are large (larger in LCY_sulf, c.f. Figure 3b), even in a five year mean, with most continental areas experiencing at least 4 K cooling, locally dropping <10 K over central Asia

Sea-ice/ocean changes for 20 years
The long lasting global cooling response cannot be explained by the direct radiative forcing from the volcanic aerosols, since the aerosols have mostly disappeared after 5 years. In Figure 3c we show the 12 month running mean change in global mean sea ice area. Sea ice immediately response to the eruption induced surface cooling with a peak increase of sea ice area globally up to 40 % in LCY_full and up to 50 % in LCY_sulf. Global sea ice area in the model experiments is not back to climatological values before at least 20 years after the eruption. When inspecting NH sea ice and ocean changes more in detail ( Figure S4) we find that Arctic sea ice area is increased immediately after the eruption and for more than 20 years with a maximum of 7 million km 2 (not shown), a 104 % increase in post-eruption year 2. This change is accompanied by a reduction of ocean heat content (not shown) and a decrease of poleward ocean heat transport at 60° N after the eruption, lasting from post-eruption year 3 up to 20 with a maximum decrease up to 0.1 PW (20 %) in post-eruption year 5 ( Figure   S4). Thus, abrupt surface cooling and decrease of upper ocean heat content in the NH leads to an immediate increase of Arctic sea ice area in the first years. The reduced poleward ocean heat transport at northern mid-latitudes for up to 20 years sustains the sea ice and climate surface cooling signal for more than 20 years in the NH and also globally. Antarctic sea ice area reveals a slightly later and shorter lasting increase from post-eruption years 1 to 5. The poleward ocean heat transport at 60°S is much more variable than in the NH and does not show significant changes over longer time periods in our simulations. This may be due to the later supply of AOD to the southern hemisphere (SH) ( Figure S1), thus later radiative forcing, as well as smaller AOD and hence weaker surface climate response in the SH compared to the NH (Figures 4, and   5). For a tropical January eruption, AOD is first distributed in the tropical belt in the first few weeks before being transported poleward to the NH winter/spring season and then to the SH in the following months ( Figure S1; see also Toohey et al., 2011), reflecting the pathways and seasonality of the Brewer Dobson circulation (Plumb, 2002). Atmospheric circulation changes are expected to be significant for the LCY eruption as was shown by Toohey et al. (2011Toohey et al. ( , 2013 and will be further investigated in a follow up paper.

Large impacts on precipitation and vegetation
A strong cooling of the atmosphere, like from an explosive volcanic eruption, leads to decreased precipitation (Robock and Liu, 1994). In our simulations, global mean precipitation (Figure 3d) decreases ~25 % (~0.8 mm/day) in the LCY_full scenario and more than ~30 % (1 mm/day) in the LCY_sulf scenario. The LCY_sulf simulation is outside the two standard deviation range of the LCY_full ensemble. Return to background climatological precipitation takes more than 10 years in both scenarios. The minimum precipitation is found between 2 and 3 years after the eruption, closely following the drop in the temperature signal. throughout the tropical East Atlantic, Northern Africa, Middle East, and the Arabian Peninsula. These are relatively dry regions, so a moderate absolute increase in precipitation (<1 mm/day) corresponds to more than a doubling of rainfall over large parts of this region. Comparing LCY_full and LCY_sulf, the impacts are generally stronger for the latter scenario both where we find drying and wetting.
The strong AOD increase, global surface cooling and decrease in precipitation together results in a decrease in land plant productivity (Net Primary Production (NPP) of >30 % during the first three years after the eruption, followed by suppressed production during the next ~15 years (Figure 3d). NPP is especially reduced over the NH land with peak decrease >75 % over high latitudes and a gradual weakening of this signal towards lower latitudes (Figures 5g,h and S3). However, over Northern Africa and surrounding areas, where precipitation increases significantly due to the southward shift of the ITCZ, we find a corresponding enhancement in land plant productivity as shown by a strong increase in NPP in this region. This is by far the strongest signal we detect in NPP with more than 400 % gain in some areas.

El Niño conditions
The ENSO response of the simulations is shown in Figure 6. Even though the initial conditions of the experimental set-up span different ENSO states, there is a rapid convergence to a robust response in the two eruption scenarios, ONI values drop below -5 during the first 5 years after the eruption. First, the LCY_full ensemble spread is suppressed for 3 years after the eruption, before beginning to diverge again. The ONI values drop far below the range of natural variability, but there are two local maxima during post-eruption years 0 and 1, with an amplitude similar to a strong El Niño in the control simulation.
Thus, we compare the ensemble mean of these temporal peaks to the ensemble mean average temperature over post-eruption years 0-5 and 1-4, respectively. The resulting spatial patterns are shown in Figure 6b and c. Compared to the dominant volcanic surface cooling pattern, these maxima represent El Niño conditions, indicating that the super volcanic eruption is triggering an El Niño response in the first two years after the January eruption, which is masked by the strong surface cooling caused by the eruption.

Comparison with other studies
Our simulations show very large climate impacts from the LCY sulfur-and halogen-rich super-eruption, and a break from the tendency towards simulating smaller climate impacts from the Toba eruption, that we noted in the Introduction. In Figure   7 we show scatter plots comparing our simulations of LCY to other simulations of super volcanic eruptions using sulfur only injections (English et al., 2013;Jones et al., 2005;Metzner et al., 2014;Robock et al., 2009;Timmreck et al., 2010).
Compared to other model studies with interactive aerosols of volcanic eruptions of magnitude Mv > 7, our simulations show very large maximum AODs and thus maximum surface climate impacts for a given sulfur injection (Figure 7 a, b). The largest climate cooling for a super eruption is achieved when using linearly scaled AOD values based on Pinatubo (Jones et al., 2005;Robock et al., 2009), but this approach is simplified, since there are several feedbacks (i.e., self-limiting, scattering, and removal of aerosols) that makes the relationship between sulfur injection, aerosols, radiative forcing, and climate highly non-linear i.e. (Bekki, 1995;Metzner et al., 2014;Pinto et al., 1989).
Limiting our comparison to model studies that use sulfur injection to generate self-consistent AOD estimates, we see that our model experiments show larger radiative impacts and larger surface cooling per injected sulfur mass to the stratosphere than those studies (English et al., 2013;Metzner et al., 2014;Timmreck et al., 2010). The differences could be caused by different aerosol microphysics (bulk vs. modal vs sectional modules), radiation, advection and depositions schemes (see discussions by English et al 2013;also Marshall et al. (2018) for a discussion), as well as atmospheric chemistry (OH, ozone, H2O) and climate/ESM model differences (coupling, resolution, clouds). Our comparison (Figure 7) is limited by the fact that the simulations were not part of a coordinated model intercomparison yet, thus the model experiments are all different related to eruption strength, date, location, and injection altitude. Volcanic aerosol climate model intercomparisons are in progress now (see (Timmreck et al., 2018;Zanchettin et al., 2016)) and should include extremely large to super-size eruptions, where the model spread is even larger (Figure 7) but observational evidence is lacking.
Even though the halogen injection efficiency for a super eruption like LCY is highly uncertain, we expect that the effects of injected halogens would be qualitatively similar independent of the magnitude of the injection as our model results for 10 % and 1 % halogen injections revealed. Injecting additional volcanic halogens into the stratosphere leads to ozone depletion (this study; (Brenna et al., 2019) and the interaction with the OH availability impacts the aerosol formation leading to smaller maximum AOD and hence weaker surface cooling. Including the volcanic release of halogens as well as sulfur should be a part of future model intercomparisons focusing on volcanic impacts on climate and ozone. Overall, the model uncertainties for super eruptions is dramatically reduced when comparing the peak climate response to a given maximum AOD (Figure 7), which shows a clear relationship, slightly non-linear, with decreasing surface temperature with increasing AOD.
When atmospheric temperatures drop after a volcanic eruption, changes in energy balance of the climate system leads to decreased global mean precipitation (Iles et al., 2013;Robock and Liu, 1994). The global mean precipitation response to the super volcanic eruptions follows a nearly linear relation with temperature (Figure 7d), larger cooling leads to larger negative precipitation anomalies through a weakening of the global hydrological cycle since lower temperatures leads to lower relative humidity in the troposphere. Our modeled southward shift of the ITCZ towards the southern tropics is accompanied by increased precipitation across North Africa and the Middle East, which is partly simulated also in Robock et al. (2009) and Timmreck et al. (2010Timmreck et al. ( , 2012, but the area experiencing wetting is larger in our simulations. The wetting of North Africa and the Middle East in our simulations leads to a strong increase in NPP in this area and thus likely to a greening of the Sahara. Timmreck et al. (2012), with only half the volcanic forcing that we simulate, shows NPP maps (vegetation impacts are simulated using an off-line vegetation model), and here there is very little change over the first three post-eruption years throughout this region. While we cannot compare our NPP field directly to the changes in vegetation types presented in Robock et al. (2009), we note that they simulate an increase in grass cover throughout North Africa and parts of the Middle East where there is very low vegetation cover in their control run, which would imply an increase in NPP in this region as well.
We simulate El Niño conditions to our LCY super eruption during the first two post-eruption years, but superposed on a strong surface cooling signal. El Niño conditions may be favored up two years as discussed in more detailed by Emile-Geay et al. (2008) due to the uniform solar dimming leading to a thermostat mechanism (Clement et al., 1996)  Atmospheric circulation changes (i.e., stratospheric polar vortices, Annular Modes) are expected to be significant as was shown by Toohey et al. (2011Toohey et al. ( , 2013 and will be further investigated in a follow up paper. Using a fully coupled ESM with interactive aerosols and atmospheric chemistry is currently the best possible way to simulate the impacts of super-volcanoes on the Earth system. Our model setup takes into account the interactive coupling between most of the components of the Earth system, including ocean, sea-ice, bio-geochemistry, land surface and vegetation interactions. In addition, the inclusion of interactive aerosols and atmospheric chemistry is crucial to correctly simulate the feedbacks between the chemical composition of the atmosphere, aerosols and radiation. That said, there is still considerable uncertainty in the impacts of volcanic sulfur injections, particularly in the conversion of SO2 into sulfate aerosol, aerosol size and the lifetime of the radiative perturbation. The uncertainty in the Earth system's reaction to a given volcanic aerosol radiative forcing seems to be smaller (c.f. Figure 7c). Based on our set of model experiments and sensitivity studies next to existing model studies on the volcanic impact on ozone (Black et al., 2014;Brenna et al., 2019;Klobas et al., 2017;Muthers et al., 2015;Robock et al., 2009), we evaluate the ozone response for volcanic sulfur and halogen injections into the stratosphere currently to be more robust than for the climate response to a given sulfur injection.
Finally, LCY may be the super eruption of the last 100 kyrs as a slightly higher sulfur mass loading was just released (Cisneros et al. to be submitted) and Toba is estimated to be less sulfur rich than previously assumed (Chesner and Luhr, 2010). To compare these two super eruptions and petrological estimates other archives such as ice core records would be needed. However, no tephra has been identified in Greenland and Antarctica ice cores for both eruptions up to now (Abbott et al., 2012;Svensson et al., 2013). This model study together with the new dating of the LCY eruption to 80.8±6.7 kyrs and a higher (Cisneros et al., to be submitted) will hopefully stimulate upcoming studies finding corresponding paleo proxies in ice cores, climate, and archaeological archives.

Summary and conclusions
Our comprehensive aerosol chemistry Earth System Model (ESM) study shows that a sulfur-and halogen-rich tropical super-eruption like Los Chocoyos (LCY) has massive impacts on global climate and the environment lasting at least 20-30 years.
In the model, enhanced volcanic sulfate burdens and aerosol optical depth (AOD) persists for five years, while the halogens stay elevated for ~14 years. Under pre-industrial conditions, the eruption leads to a global collapse of the ozone layer (80 % decrease) with global mean values of 50 DU and increasing surface UV-B by 550 % globally over the first 5 years after the eruption. (In high latitudes the increase is >1000 %). The ozone layer takes 15 years to recover. The simulated volcanic eruption, at 15° N in January, shows an asymmetric hemispheric response with enhanced AOD, ozone, UV, and climate signals over the Northern Hemisphere (NH).
The eruption cools the global climate lasting more than 30 years with the peak AOD of >6 leading to surface cooling >6 K and precipitation and terrestrial net primary production (NPP) decreases up to 30 % in the first two years. Locally, a wetting (>100 %) and strong increase of NPP (>400 %) over Northern Africa is simulated in the first five years related to a southwards shift of the Inter-Tropical Convergence Zone to the southern tropics. Global sea ice area almost doubles, and the long lasting surface cooling is sustained by sea ice/ocean changes mainly in the Arctic showing an immediate sea ice area increase followed by a decrease of poleward ocean heat transport at 60° N from year 3, both changes lasting up to 20 years.
The ocean responds with El-Niño conditions in the first two years which are masked by the strong volcanic induced surface cooling.
The long lasting surface cooling is sustained by sea ice/ocean changes showing an immediate sea ice area increase followed by a decrease of poleward ocean heat transport at 60° latitude from year 3, both changes lasting up to 20 years.
In contrast, simulations of LCY including sulfur, but neglecting halogens, reveal larger sulfate burden and maximum AOD (~8), hence a larger radiative forcing with more pronounced surface climate cooling (>7 K) and reduced precipitation (25 %) globally, even though spatial patterns of changes are similar to the simulations including volcanic sulfur and halogens. The environmental impacts reveal the opposite signal with a short-lived increase of column ozone of 100 DU (>30 %) and decrease of UV (>60 %) lasting up to 3 years.
LCY is one of the largest volcanic eruptions over the past 100 kyrs and we predict large impacts on the biosphere and thus any human populations at that time. Finding paleo proxies showing the impact of LCY on climate and the environment should be possible, given the large long lasting impact from our ESM simulations, but will require high (decadal) temporalresolution archives using the eruption as a time marker.

Code and data availability
All simulation data will be archived in the NIRD Research Data Archive on acceptance of the manuscript. Post processing and visualization of data was performed with Python and the code and post processed data files are available on request from the corresponding author.

Supplement
The supplement related to this article is available online at:

Author contributions
HB performed the simulations, data analysis and produced the figures. HB, KK and SK interpreted the results. MM provided the CESM2(WACCM) model and supported HB in performing simulations. HB wrote the manuscript with contributions from all co-authors.