Decadal Disruption of the QBO by Tropical Volcanic Supereruptions

The Los Chocoyos (14.6°N, 91.2°W) supereruption happened ∼75,000 years ago in Guatemala and was one of the largest eruptions of the past 100,000 years. It emitted enormous amounts of sulfur, chlorine, and bromine, with multi‐decadal consequences for the global climate and environment. Here, we simulate the impact of a Los Chocoyos‐like eruption on the quasi‐biennial oscillation (QBO), an oscillation of zonal winds in the tropical stratosphere, with a comprehensive aerosol chemistry Earth System Model. We find a ∼10‐year disruption of the QBO starting 4 months post eruption, with anomalous easterly winds lasting ∼5 years, followed by westerlies, before returning to QBO conditions with a slightly prolonged periodicity. Volcanic aerosol heating and ozone depletion cooling leads to the QBO disruption and anomalous wind regimes through radiative changes and wave‐mean flow interactions. Different model ensembles, volcanic forcing scenarios and results of a second model back up the robustness of our results.

Extremely large volcanic eruptions which inject several 100 Tg SO 2 might therefore have the potential to disturb the QBO and may shed light on the question of whether the volcanic response of the QBO is indeed different from the geoengineering scenarios.
Recently, the QBO was disrupted in 2015-2016 (Newman et al., 2016;Osprey et al., 2016) and again in 2019-2020. Several studies investigated the possible reasons for the anomalous QBO behavior in 2015-2016, which started with reversed westerly zonal wind at 40 hPa (e.g., Barton & McCormack, 2017;Coy et al., 2017;Li et al., 2020) and similarly in 2019-2020 (Anstey et al., 2020, preprint). Their results showed that Rossby waves propagating from the extratropical Northern and Southern Hemisphere have been the most likely cause of the unusual behavior.
Here, we present a set of sulfur-and halogen-rich Los Chocoyos-like supereruption experiments with two different aerosol climate models, both including an interactive QBO in the tropical stratosphere. The climate and environmental impacts of the Los Chocoyos simulations were previously discussed in Brenna et al. (2020). In this paper, we study the impacts of volcanic eruptions on the stratospheric circulation and especially focus on the potential dynamics of supereruptions in modifying the QBO, how long the impacts last, and how the tropical circulation may be influenced.

CESM2(WACCM6)
In this study, we run the fully coupled Community Earth System Model, Version 2 (CESM2) (Danabasoglu et al., 2020). The atmospheric component is the Whole Atmosphere Community Climate Model (WACCM6), based on the Community Atmospheric Model Version 6 (Gettelman et al., 2019). We refer to this setup as CESM2(WACCM6). A comprehensive ESM with interactive aerosols and atmospheric chemistry spanning the whole atmosphere from the surface to the lower thermosphere with model top at ∼140 km altitude. The middle atmosphere chemistry module includes the SO x , O x , NO x , HO x , ClO x , and BrO x chemical families, implementing 98 compounds and ∼300 reactions including gas-phase, photolytic and heterogeneous reactions on aerosols. Stratospheric sulfate aerosols are formed interactively from sulfur compounds and modeled by the modal aerosol model MAM4 (Liu et al., 2016;Mills et al., 2016). The atmospheric horizontal resolution is 0.95° longitude by 1.25° latitude with 70 hybrid sigma pressure layers. CESM2(WACCM6) has parameterizations for orographic (Scinocca & McFarlane, 2000) and non-orographic (both frontal and convective; Richter et al., 2010) gravity waves. The QBO is internally generated.
The ocean model of CESM2(WACCM6) is the Parallel Ocean Project v. 2 (POP2) model running at ∼1 × 1° horizontal resolution with 60 layers in the vertical. The sea ice model is CICE5 (Bailey et al., 2018), which runs on the same horizontal grid as POP2. The land surface model is the community land model Version 5 (CLM5), which runs on the same horizontal grid as WACCM6 with dynamic vegetation, interactive biogeochemistry (carbon, nitrogen, methane), and prognostic crops (Fisher et al., 2019). We run CESM2 (WAC-CM6) under 1850 pre-industrial conditions.

MAECHAM5-HAM
The aerosol microphysical model HAM (Stier et al., 2005) is interactively coupled to the general circulation model (GCM) MAECHAM5 (Giorgetta et al., 2006) and was extended to a stratospheric version (Niemeier et al., 2009). The model runs with a horizontal spectral resolution of T42 and 90 vertical layers up to 0.01 hPa (~80 km altitude) and internally generates a QBO (Niemeier et al., 2020). The parameterization of the momentum flux deposition for non-orographic atmospheric gravity waves as used in MAECHAM5-HAM follows Hines (1997aHines ( , 1997b. The parameterization was changed to space-and time-constant gravity wave parameters for the source spectrum (Manzini & McFarlane, 1998;Roeckner et al., 2003). MAECHAM5-HAM includes a simplified sulfur chemistry to simulate SO 2 oxidation, but no complex chemistry for the precursors. Therefore, OH, ozone, NO x etc., concentrations are prescribed on a monthly mean basis (Timmreck et al., 2003). The limitation of OH in case of a supereruption is parameterized after Bekki (1995); see Timmreck et al. (2010) for details. The model is set up for present day conditions with monthly mean climatological sea surface temperatures.

Model Experiments
To model the impact of the Los Chocoyos-like eruption on the atmosphere, we use the published petrologically estimated erupted sulfur and halogen masses (523 Mt S, 1,200 Mt Cl, 2 Mt Br) from Brenna et al. (2020) as input. We inject all of the erupted sulfur mass and 10% of the erupted chlorine and bromine mass into the stratosphere, a conservative estimate (see the detailed discussion in Kutterolf et al., 2013, Krüger et al., 2015and Brenna et al., 2019. SO 2 , HCl, and HBr are injected into the model grid boxes at 14.6°N at 24 km altitude. This injection level was chosen as a best estimate following the observed spread of the volcanic cloud after the Pinatubo eruption (Guo et al., 2004). To technically enable this supereruption model run, we had to spread the emissions of volcanic gases over adjacent model columns and time steps. We chose to spread between 80° and 97.5°W over the six first days of January, which is a reasonable estimate for the minimum duration of the eruption based on the given erupted mass of the eruption and the respective magma discharge rate (Rose et al., 1987;Wilson & Walker, 1987).
The model experiments are the same as in Brenna et al. (2020) and are summarized in Table 1. We run an ensemble of six simulations starting from different ENSO and QBO states of the preindustrial control run with the combined sulfur and halogen forcings (LCY_full). In addition, we perform two simulations (ENSO neutral: QBO easterly and westerly) with only the sulfur forcing (LCY_sulf) to explore the response between the different volcanic forcing scenarios. All ensemble members last 35 years, and the control simulation is 70 years with constant pre-industrial (1850) forcings. We have also 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 (LCY_1%halog) instead of 10%. A single simulation using MAECHAM5-HAM corresponding in forcing to the LCY_sulf.qboW simulations is used to compare the QBO response across the two models.

QBO in the Control Simulation
The QBO in CESM2(WACCM6) has previously been evaluated in Gettelman et al. (2019). The QBO in our control simulation has a period of approximately 27 months at 30 hPa varying between 20 and 33 months ( Table 1, Figures 1a and 1b). The modeled QBO winds are in thermal wind balance with the tropical temperature field, as shown by the wind shear (supporting information Text S1 and Figures 2a and 2b). We also see that the alternating patterns of easterly and westerly zonal wind shears align with the alternating patterns of positive and negative monthly mean temperature anomalies (Figures 1a and 1b; all anomalies are with respect to the monthly climatology). The thermal wind balance relates the shear of the zonal wind to the meridional temperature gradient and is a valid approximation under geostrophic balance (Baldwin et al., 2001).
The forcing terms of the QBO can be diagnosed using the Transformed Eulerian Mean (TEM) zonal momentum equation (Andrews & McIntyre, 1976, see supporting information Text S2). Consistent with the QBO theory (Baldwin et al., 2001), we find that resolved large scale waves such as Kelvin and Rossby-gravity waves and unresolved waves such as parameterized small scale orographic and non-orographic gravity waves deposit their momentum in the wind shear zones (Figure S1i-1j), which causes the winds to descend.

Super Volcanic Eruption Disrupts the QBO
The sulfur and halogen injection from the Los Chocoyos-like supereruption disrupts the QBO for ∼10 years (Figures 1c-1d). The equatorial winds in the lower stratosphere are forced into prolonged easterly winds for about 5 years, independent of the initial QBO and ENSO phases, followed by prolonged westerlies. The tropical stratospheric circulation resumes to a QBO regime after a decade.
The injection of sulfur and halogens to the stratosphere leads to a very strong warming in the tropical lower and middle stratosphere caused by increased aerosols in this region (Figures 3b,. This warming lasts approximately 4 years. In the middle and upper stratosphere, volcanically induced ozone depletion (Figures 3a and S2a) leads to cooling which intensifies over the first five post-eruption years (Figures 1c  and 1d). The cooling lasts in total up to 10 years in the tropics until halogens and ozone have recovered to background levels (see also Brenna et al., 2020).
Even though the disruption to the QBO system is very strong, the thermal wind balance still holds to the same extent as in the control simulation (Figures 2a-2d). The shear in the volcanically forced easterlies and westerlies correspond directly to changes in the meridional temperature gradient in the tropics. The vertical temperature anomaly structure and the zonal wind shear are decoupled with strong cooling coinciding with the region of change from easterly to westerly shear during post-eruption years 3-6. This holds for both the initial easterly and westerly ensembles (Figures 2c-2d). When the tropical winds shift to prolonged westerlies (years 4-9), a similar response with reversed signals is observed.
After the resumption of the QBO regime around year 10, the period of the oscillation is still different from the CTR. Using Fourier analysis on the last 25 years of each simulation in the LCY_full ensemble, we find a QBO period between 30 and 33 months in each simulation (Table 1), compared to 27 in the CTR when analyzing 25-year slices. This shows that the impact on the QBO persists longer than the volcanic forcing lasts. This may be related to the strong surface cooling and reduced hydrological cycle in the tropics which lasts more than 30 and 20 years, respectively (Brenna et al., 2020). Inspecting the vertical flux of zonal momentum from resolved tropical waves (taken from Eq. SE5; not shown here) reveals, however, no clear change in the magnitude through the tropopause (at 100 hPa) after 10 years.

QBO Response to Volcanic Forcing
To analyze the QBO response to the volcanic forcing, we inspect the temporal-height evolution of a single ensemble member from the LCY_full in more detail in Figure 3 (LCY_full.qboW.ensoN in Table 1). This single ensemble member is representative of the ensemble mean as the model spread between the ensemble members is very low for the first 10 years (not shown). For comparison, we display the same quantities for the CTR and LCY_full.qboE.ensoN in Figures S1-S2. Figure 4 shows the meridional-height plane of the LCY_full ensemble mean averaged for the volcanic induced easterly and westerly periods.
The large positive aerosol surface area density anomaly (Figure 3a) lasts for 5-6 years and is the dominating volcanic forcing leading to radiative warming in the lower to middle stratosphere. Ozone loss is evident throughout the whole stratosphere ( Figure 3b) and lasts for more than 10 years due to the presence of vol-BRENNA ET AL.
10.1029/2020GL089687 5 of 13 canic halogens (Brenna et al., 2020). The ozone loss peaks in the middle and upper stratosphere associated with radiative cooling.
The diabatic heating (Figure 3c), which is dominated by radiative effects (not shown here), is balanced in the TEM thermodynamic equation (Eq. SE2) by the vertical advection of heat ( Figure 3d) through a thermal direct cell. These two terms together create the important signature of the temperature tendency evolution in the simulation (Figures 3e and 3f). Through the thermal wind balance, the meridional temperature gradient controls the vertical wind shear, and the acceleration terms of the TEM momentum equation interact with the temperature field by changing the vertical wind shear. During the volcanically induced easterly phase, the meridional temperature gradient in the tropics is weak in both hemispheres, positive on the NH BRENNA ET AL.
10.1029/2020GL089687 6 of 13 and negative on the SH (not shown here) leading to weak easterly shear during the first four years from the middle to the lower stratosphere (gray contours in Figures 3 and 2d).
Gravity and resolved equatorial wave activity are impacted, as apparent in the anomalous easterly acceleration of equatorial gravity wave drag and EP flux divergence from the upper troposphere to the middle stratosphere (∼200-5 hPa) during the first 3 years following the LCY eruption (compare Figures 3 and S2 with CTR in Figure S1). The parameterized gravity wave drag contributes with weak easterly forcing on the zonal wind between 100 and 7 hPa and a very strong westerly forcing above this layer during the first two post-eruption years (Figure 3i). This westerly acceleration propagates downward together with the westerly shear zone in a step-like fashion, caused by wave activity and thermal forcing changes acting together, contributing to the zonal wind change after ∼4-5 years. Inspecting the vertical flux of zonal momentum from resolved tropical waves through the tropopause reveals a reduction (not shown here) related to reduced convection due to surface or tropospheric cooling. Further analyses of wave frequencies and characteristics combined with the lowering of the tropical tropopause height (Figures 1-3 and S2) and the shift of the ITCZ to the southern tropics (Brenna et al., 2020) would be interesting to study in future work.
The contribution by the resolved waves (mainly large-scale Kelvin and Rossby-gravity waves) is mostly in phase with the gravity wave drag below 10 hPa and opposing above (Figure 3j). Since the waves deposit their momentum in the shear zones, the very weak shear, together with the reduced wave momentum flux into the stratosphere during the first two post-eruption years (Figure 2), contributes to the weak wave forcing between 100 and 10 hPa in the same period (Figures 3i and 3j). The main opposing acceleration forcing is the vertical advection of easterly zonal momentum (Figure 3h), which tends to prevent the downward propagation of the westerly shear zone. Below the level of strong gravity wave forcing (between 7 and 1 hPa), a strengthened secondary meridional circulation results with more upward and northward residual flow (Figures 3g, 4i and 4k).
Comparing temperature and velocity tendencies calculated using the TEM equations (Figures 3e and 3k) with tendencies calculated directly from the model temperature and velocity fields (Figure 3f and 3l) reveals that most of the evolution of the modeled tendency fields can be explained using the theoretical TEM framework (supporting information Text S2).
During post-eruption years 5-9, the aerosol forcing has mostly disappeared and the ozone forcing is decaying. Changes in the meridional temperature structure in the lower to middle stratosphere leads, through the thermal wind balance, to a westerly phase of the zonal wind shear from the middle to the lower stratosphere. A new easterly phase is forming above and is descending in the same step-wise fashion forced by the wave breaking in the shear zone.
The volcanically induced wind changes in LCY_full are accompanied by changes in the residual circulation with first an increased meridional overturning circulation in the tropical stratosphere toward the Northern Hemisphere during the easterly wind period, which is weakened in the following westerly wind period (Figure 4). In the LCY_sulf ensemble, which is forced only by sulfur injection, we see the same transition to volcanic easterly zonal winds at the equatorial stratosphere (Figures 1e-1f). The duration of this induced easterly wind is approximately the same as in LCY_full. This is followed by a westerly wind phase, which is not as prolonged and strong as in LCY_full, lasting only 3 years.
The relationships between the TEM thermodynamic and momentum forcing terms, the zonal wind shear and the temperature and momentum tendencies are similar in LCY_full (Figures 3 and 4, S2) to LCY_sulf ( Figures S3-S5). The major difference is the signal and length of the forcing from the volcanically induced ozone response. The ozone loss associated cooling in LCY_full lasts more than 10 years in the tropics (Figure 4), whereas in LCY_sulf the main signal is ozone increase maximizing in the upper stratosphere with a small ozone decrease in the mid stratosphere ( Figure S5). There is only a relatively small and relatively constant (in time) ozone increase induced warming perturbation left between post-eruption years 5-8.
BRENNA ET AL.

Comparison with MAECHAM5-HAM and Other Model Sensitivity Tests
To check whether the QBO response during the first 1-2 years is model dependent, we compare results for the sulfur only and QBO W experiments with results of MAECHAM5-HAM using a similar experimental set up ( Figure S6; see Methods). In MAECHAM5-HAM we also simulate an anomalous interruption of the westerly wind and shift to an easterly phase connected to the volcanic aerosol forcing in the lower stratosphere. The easterly wind regime is prolonged and last almost 2 years. This contrasts with the ∼5 years length of easterly wind regime in WACCM6 and is related to the much shorter aerosol forcing in MAECHAM5-HAM, maximum within 4-9 months, turning into a cooling tendency in the lower stratosphere afterward. Most of the volcanic induced temperature anomaly has dissipated after 1.5 years.
Previous comparisons of WACCM and ECHAM for volcanic (Brenna et al., 2020;Marshall et al., 2018;Zanchettin et al., 2016) and geo-engineering (Niemeier et al., 2020) studies revealed a similar picture: MAECHAM5-HAM generates peak AOD more rapidly and also decays faster than CESM2(WACCM6). This is at least partly related to interactive atmospheric chemistry taken into account in CESM2(WAC-CM6), where OH limitation delays the aerosol formation and leads to significantly smaller aerosols in WACCM than in MAECHAM5-HAM for super-size eruptions (Brenna et al., 2020;Timmreck et al., 2010). Even though the aerosols behave very differently in the two models, the peak tropical AOD and maximum temperature anomalies in the equatorial belt are of similar magnitudes leading to the same easterly zonal wind response in the tropical stratosphere due to a sulfur rich super volcanic eruption ( Figure S6).
The pattern of temperature anomalies and the abrupt transition from westerly to easterly wind is very similar between the two model simulations if the different internal time scales for the aerosol forcing in the two models is taken into account. On the other hand, the long-term responses in WACCM and ECHAM are diverging due to a different volcanic forcing duration next to the lack of an interactive ocean model in the ECHAM model.
Additionally, we also tested the sensitivity to different volcanic forcing scenarios with CESM2(WACCM6). Decreasing the halogen injection efficiency to 1% in LCY_1%halog compared to the 10% of LCY_full, but keeping the sulfur forcing the same, results in a volcanic QBO response which is intermediate between the LCY_full and LCY_sulf responses (Figure 1g). The volcanic induced easterly winds are very similar to the one in LCY_full, while the following westerly phase is shorter and more similar to LCY_sulf. Decreasing the SO 2 injection to 10 Mt results in a much weaker QBO impact, as seen in LCY_1%sulf (Figure 1h). Similar to what was observed after Pinatubo (Labitzke, 1994), injecting between 10 and 20 Mt SO 2 (Timmreck et al., 2018) led to a delay in the descent of the following easterly phase.
Even though the temperature anomaly patterns between the LCY_full and LCY_sulf ensembles are very different, the disruption of the QBO is very similar during the first post-eruption years. This indicates that the initial (change to easterly phase) QBO disruption is mainly driven by the volcanic aerosol heating in the lower to middle stratosphere in the core region of the QBO. The temperature tendency is quite similar in this region between the two simulations indicating that it is a more important driver of the disruption than the temperature anomaly.
The threshold for triggering a disturbance of the QBO by a ∼15°N volcanic eruption in CESM2(WACCM6) is around the size of the Pinatubo eruption. As noted above, we find a small disturbance in the LCY_1%sulf injection, which is consistent with the data after Pinatubo (Labitzke, 1994). However, the threshold value might be highly model-dependent as a multi-model geoengineering study has revealed (Niemeier et al., 2020). Further model simulations, including model intercomparisons are needed to clearly identify the threshold value.
Overall, super volcanic eruptions in the tropics, either rich in sulfur and halogens or in sulfur only, disrupt the QBO and lead to anomalous easterly winds in the tropical stratosphere during the length of the volcanic forcing based on different model ensemble members and different interactive aerosol climate models. This result stands in sharp contrast to geo-engineering model experiments which reveal a continuously westerly wind regime response for artificial SO 2 injections into the equatorial stratosphere (Aquila et al., 2014;Niemeier & Schmidt, 2017;Richter et al., 2017). In contrast to our eruption simulations, the geo-engineering experiments continuously injected aerosols or source gases at the equator, producing a smaller but continuous aerosol forcing and a different meridional temperature gradient with the highest temperatures centered on the equator. Consequently, the impact on the thermal wind balance differs and hence westerly winds are developing in the corresponding atmospheric levels. Richter et al. (2017) showed that this disruption is not simulated in CESM1(WACCM) for geo-engineering experiments where injections occur at 15° latitude on both hemispheres, which is closer to our eruption experiments than the equator.
This study is just a first step into this research topic. More detailed investigations are needed in future studies, for example, taking different volcanic eruption locations in the tropics into account. Furthermore, the radiative effects of SO 2 may also play significant roles (Osipov et al., 2020).
Our model results are not only relevant for the distribution and deposition of volcanic ash and gases of the LCY or a similar eruption in the tropics but are also important for a better understanding of the QBO phenomenon in itself. Recent observations revealed an unusual behavior of the QBO in 2015-2016 with an anomalous upward displacement of the westerly phase and easterly winds below, which were related to the easterly momentum of extratropical planetary waves in the tropics (Newman et al., 2016;Osprey et al., 2016) and similar in 2019-2020 again (Anstey et al, 2020). In contrast to these events, our QBO disruption starts as a result of intense volcanic aerosol heating and subsequent dynamical changes in the wave mean flow interaction in the QBO core region in the tropics. Future sensitivity tests and model intercomparisons (e.g., Butchart et al., 2018) will be useful to further investigate QBO processes and to explore the sensitivity of the QBO response to different volcanic forcing.

Summary and Conclusions
In this paper, we have simulated the impact of a sulfur-and halogen-rich supereruption like the one of Los Chocoyos (14.6°N, 91.2°W) on the QBO with the Earth System Model CESM2(WACCM6). The model includes aerosol formation, interactive atmospheric chemistry, and an internally generated QBO.
The model results reveal that such a supereruption disrupts the QBO for roughly a decade ( Figure S7). The equatorial winds in the lower to middle stratosphere are forced into prolonged easterly winds, independent of initial conditions, followed by prolonged westerlies. The external forcing from aerosols (heating) and aerosol induced ozone depletion (cooling) in a halogen rich atmosphere leads to changes in the temperature of the tropical stratosphere, impacting the zonal wind shear through the thermal wind balance. This changes the wave propagation and suppresses the downward propagation of the westerly phase in our model set up for the first five post eruption years. During the following 5 years, anomalous westerly winds develop. This change is driven by thermodynamical changes due to the decaying volcanic aerosol and ozone depletion forcing in the middle to upper stratosphere in combination with wave forcing. Following the major QBO disruptions of the first 10 post-eruption years, the volcanic forcing in the stratosphere decays significantly, and QBO conditions return but with a markedly longer period (prolonged from 27 to more than 30 months).
We tested the robustness of our model response to a supereruption at ∼15°N by running several ensemble members with different QBO and ENSO starting conditions and different volcanic forcing scenarios which gave similar results. The abrupt onset of anomalous easterly winds is a robust feature of the same experimental set-ups in the WACCM6 and the MAECHAM5-HAM models with sulfur only injections. The signal, pattern, and mechanism of changes are similar, but the duration of the volcanic forced easterlies differs. This is related to the different lifetimes of stratospheric aerosols between the two models and partly due to the role of interactive atmospheric chemistry.
Overall, we conclude that super volcanic eruptions in the tropics in the past, either rich in sulfur and halogens or in sulfur only, can disturb the QBO in the tropical stratosphere by initiating easterly winds after the first months of the eruption with a corresponding enhanced secondary meridional circulation. The presence of volcanic chlorine and bromine leads to ozone depletion and thus cooling, whereas a sulfur-only paleo-eruption leads to ozone increases and thus heating in the middle to upper stratosphere.
Our results are relevant for the transport of volcanic material after a supereruption in the tropics, the interpretation of volcanic supereruption markers in paleo proxy records, and may lead to a better understanding of the QBO.

Data Availability Statement
CESM2 is open software and available from NCAR. All CESM2 simulation data have been archived in the NIRD Research Archive (DOI: 10.11582/2020.00043) Python scripts and data required to reproduce our results are available through GitHub/Zenodo (https://doi.org/10.5281/zenodo.4028819). providing the CESM2(WACCM6) model code and for their technical model support. The CESM project is supported primarily by the National Science Foundation (NSF). This material is based upon work supported by NCAR, which is a major facility sponsored by the NSF under Cooperative Agreement No. 1852977. The CESM2 simulations for this study were performed on resources provided by UNINETT Sigma2 -the National Infrastructure for High Performance Computing and Data Storage in Norway. The MAECHAM5-HAM simulations were performed on the computer of the Deutsches Klima Rechenzentrum (DKRZ). HB was partly financed by the NFR project KeyCLIM (project number 295046). UN and CT are financed by the DFG Research Unit VollImpact (FOR2820, DFG No.398006378) within the projects VolARC (UN) and VolClim (CT).