Leon Hermanson

Robert Bilbao

Nick Dunstone

Martin Ménég

Pablo Ortega

Holger Pohlmann

Jon I Robson

Doug M Smith

Gar Stran

Claudia Timmreck

Steve Yeager

Gokhan Danabasoglu

Robust Multiyear Climate Impacts of Volcanic Eruptions in Decadal Prediction Systems

Major tropical volcanic eruptions have a large impact on climate, but there have only been three major eruptions during the recent relatively well-observed period. Models are therefore an important tool to understand and predict the impacts of an eruption. This study uses five state-of-the-art decadal prediction systems that have been initialized with the observed state before volcanic aerosols are introduced. The impact of the volcanic aerosols is found by subtracting the results of a reference experiment where the volcanic aerosols are omitted. We look for the robust impact across models and volcanoes by combining all the experiments, which helps reveal a signal even if it is weak in the models. The models used in this study simulate realistic levels of warming in the stratosphere, but zonal winds are weaker than the observations. As a consequence, models can produce a pattern similar to the North Atlantic Oscillation in the first winter following the eruption, but the response and impact on surface temperatures are weaker than in observations. Reproducing the pattern, but not the amplitude, may be related to a known model error. There are also impacts in the Pacific and Atlantic Oceans. This work contributes toward improving the interpretation of decadal predictions in the case of a future large tropical volcanic eruption.

Introduction

Volcanoes have global impacts on climate from monthly to centennial timescales [START_REF] Ding | Ocean response to volcanic eruptions in Coupled Model Intercomparison Project 5 simulations[END_REF][START_REF] Robock | Volcanic eruptions and climate[END_REF][START_REF] Timmreck | Modeling the climatic effects of large explosive volcanic eruptions[END_REF][START_REF] Timmreck | The impact of stratospheric volcanic aerosol on decadal-scale climate predictions[END_REF][START_REF] Zhong | Centennial-scale climate change from decadally-paced explosive volcanism: A coupled sea ice-ocean mechanism[END_REF]. Volcanic impacts are needed to fully explain the climate variability of the twentieth century [START_REF] Smith | Role of volcanic and anthropogenic aerosols in the recent global surface warming slowdown[END_REF]. Understanding these impacts can help with prediction and planning in the case of a volcanic eruption today. In the modern observational period, there have been three major tropical eruptions: Agung (1962), El Chichón (1982), and Pinatubo (1991).

When a volcanic eruption is powerful enough, such as the ones mentioned above, the volcanic plume reaches up into the stratosphere and aerosols formed there can remain for months to years. When the eruptions are in the tropics, the aerosols are quickly spread across all longitudes by the strong winds in the stratosphere. Once in the stratosphere, the aerosols cause local warming. This is caused by shortwave absorption and also by blocking of upward longwave radiation [START_REF] Hansen | Radiative forcing and climate response[END_REF][START_REF] Lacis | Climate forcing by stratospheric aerosols[END_REF][START_REF] Swingedouw | Impact of explosive volcanic eruptions on the main climate variability modes[END_REF]. In winter, the stratospheric warming modifies the midlatitude to polar temperature gradient, which strengthens the zonal winds at high latitudes, and this decreases the propagation of planetary waves poleward, which further increases the zonal winds [START_REF] Stenchikov | Arctic Oscillation response to the 1991 Mount Pinatubo eruption: Effects of volcanic aerosols and ozone depletion[END_REF]. The result is an increased chance of a positive North Atlantic Oscillation (NAO) after the volcanic eruption, as discussed by many authors including [START_REF] Robock | Winter warming from large volcanic eruptions[END_REF], [START_REF] Robock | Volcanic eruptions and climate[END_REF], [START_REF] Stenchikov | Arctic Oscillation response to the 1991 Mount Pinatubo eruption: Effects of volcanic aerosols and ozone depletion[END_REF], [START_REF] Shindell | Dynamic winter climate response to large tropical volcanic eruptions since 1600[END_REF], and [START_REF] Stenchikov | Arctic Oscillation response to volcanic eruptions in the IPCC AR4 climate models[END_REF]. Note that there are many other influences on the NAO, such as the tropical Pacific, the North Atlantic, the quasi-biennial oscillation, and solar variability [START_REF] Dunstone | Skilful predictions of the winter North Atlantic Oscillation one year ahead[END_REF]. [START_REF] Ding | Ocean response to volcanic eruptions in Coupled Model Intercomparison Project 5 simulations[END_REF] found only a weak relationship between volcanoes and the El Niño-Southern Oscillation (ENSO), but more recent studies with newer models consistently find an increased chance of an El Niño within 2 years after the eruption [START_REF] Khodri | Tropical explosive volcanic eruptions can trigger El Niño by cooling tropical Africa[END_REF][START_REF] Predybaylo | Impacts of a Pinatubo-size volcanic eruption on ENSO[END_REF][START_REF] Stevenson | El Niño like" hydroclimate responses to last millennium volcanic eruptions[END_REF]. There 10.1029/2019JD031739 is also compelling evidence of this from paleoreconstructions (Adams et al., 2003;[START_REF] Wilson | Reconstructing ENSO: The influence of method, proxy data, climate forcing and teleconnections[END_REF]. One proposed mechanism is related to differential cooling between land and ocean in the West Pacific and associated changes in the Walker circulation, which are amplified by Bjerknes feedback [START_REF] Swingedouw | Impact of explosive volcanic eruptions on the main climate variability modes[END_REF].

There are also regional impacts for North Atlantic Ocean dynamics. Over decadal timescales, the Atlantic Meridional Overturning Circulation (AMOC) increases in strength and leads to regional warming [START_REF] Iwi | Mechanisms linking volcanic aerosols to the Atlantic Meridional Overturning Circulation[END_REF][START_REF] Ottera | External forcing as a metronome for Atlantic multidecadal variability[END_REF][START_REF] Stenchikov | Volcanic signals in oceans[END_REF][START_REF] Swingedouw | Impact of explosive volcanic eruptions on the main climate variability modes[END_REF]. Decreased temperature and increased salinity lead to an enhanced AMOC, and the strength of this enhancement is correlated with background variability of overturning in the model [START_REF] Ding | Ocean response to volcanic eruptions in Coupled Model Intercomparison Project 5 simulations[END_REF], though changes in windstress curl are also likely to play a role [START_REF] Mignot | Volcanic impact on the Atlantic Ocean over the last millennium[END_REF]. In some models, background conditions can alter the response in the North Atlantic [START_REF] Swingedouw | Bidecadal North Atlantic ocean circulation variability controlled by timing of volcanic eruptions[END_REF][START_REF] Zanchettin | Background conditions influence the decadal climate response to strong volcanic eruptions[END_REF].

There has been discussion on whether Coupled Model Intercomparison Project phase 5 (CMIP5) climate models are able to capture the dynamic response to volcanic aerosols (e.g., [START_REF] Driscoll | Coupled Model Intercomparison Project 5 (CMIP5) simulations of climate following volcanic eruptions[END_REF]. However, experiment design is important, the Pacific response may depend on capturing the phase of ENSO at the time of the eruption [START_REF] Lehner | The importance of ENSO phase during volcanic eruptions for detection and attribution[END_REF]. Commonly used seasonal averages may not capture the response [START_REF] Barnes | Robust wind and precipitation responses to the Mount Pinatubo eruption, as simulated in the CMIP5 models[END_REF][START_REF] Predybaylo | Impacts of a Pinatubo-size volcanic eruption on ENSO[END_REF]. Including an eruption that the model responds weakly to will reduce the signal in a multivolcano ensemble (Bittner et al., 2016). Most importantly, model experiments need to be based on large ensemble experiments to detect a small signal with respect to large internal variability (Ménégoz et al., 2018). The problem may be further exacerbated by models having a smaller signal than the real world, especially in the North Atlantic, so very large ensembles (more than a 100 members) are likely necessary to capture a response [START_REF] Baker | An intercomparison of skill and overconfidence/underconfidence of the wintertime North Atlantic Oscillation in multimodel seasonal forecasts[END_REF][START_REF] Dunstone | Skilful predictions of the winter North Atlantic Oscillation one year ahead[END_REF][START_REF] Scaife | A signal-to-noise paradox in climate science[END_REF].

Volcanic aerosols are included in decadal climate prediction and increase the skill of these predictions, particularly over the oceans (where cooling has long persistence) and Europe (Ménégoz et al., 2018;[START_REF] Timmreck | The impact of stratospheric volcanic aerosol on decadal-scale climate predictions[END_REF]. A prediction system is initialized and therefore can be used to understand better how differences in the Atlantic Ocean (Ménégoz et al., 2018) and Pacific Ocean [START_REF] Illing | Assessing the impact of a future volcanic eruption on decadal predictions[END_REF] can alter the impacts of an eruption.

In this paper, we examine the robust response, which is common to multiple models and volcanic eruptions, and (where possible) compare this with observations. Many multimodel studies use historical integrations to find the impact of volcanoes (Bittner et al., 2016;[START_REF] Driscoll | Coupled Model Intercomparison Project 5 (CMIP5) simulations of climate following volcanic eruptions[END_REF][START_REF] Iles | The global precipitation response to volcanic eruptions in the CMIP5 models[END_REF][START_REF] Khodri | Tropical explosive volcanic eruptions can trigger El Niño by cooling tropical Africa[END_REF][START_REF] Maher | Effects of volcanism on tropical variability[END_REF][START_REF] Zambri | Winter warming and summer monsoon reduction after volcanic eruptions in Coupled Model Intercomparison Project 5 (CMIP5) simulations[END_REF]. Here, we use experiments that have been run specifically to show the impact of volcanic aerosols in the context of climate predictions and are designed to reduce model biases. The models and experiments are described in the next section. In the analysis that follows, we consider the surface and zonal-mean responses first, after which we focus on the large climate modes, NAO, ENSO, and AMOC, as other impacts can often be derived from them [START_REF] Clark | Skilful seasonal predictions for the European energy industry[END_REF].

Data and Methods

This study uses decadal climate predictions that employ coupled atmosphere-ocean-ice models, which are initialized from observations and have representations of all the main external forcings. Five models are included, obtained directly from the different modeling centers and chosen because they had completed the experiments described below at the time of gathering the data for this study. The models are the European Community Earth system model (EC-Earth) V2.3 [START_REF] Hazeleger | EC-Earth V2.2: Description and validation of a new seamless earth system prediction model[END_REF]Ménégoz et al., 2018), the Max Planck Institute Earth System Model (MPI-ESM, [START_REF] Giorgetta | Climate and carbon cycle changes from 1850 to 2100 in MPI-ESM simulations for the Coupled Model Intercomparison Project phase 5[END_REF][START_REF] Pohlmann | Improved forecast skill in the tropics in the new MiKlip decadal climate predictions[END_REF], the Hadley Centre Global Environmental Model Global Coupled Configuration 2 (HadGEM3-GC2, [START_REF] Dunstone | Skilful predictions of the winter North Atlantic Oscillation one year ahead[END_REF][START_REF] Williams | The Met Office Global Coupled Model 2.0 (GC2) configuration[END_REF], the Community Earth System Model (CESM) 1.1 [START_REF] Hurrell | The Community Earth System Model: A framework for collaborative research[END_REF][START_REF] Yeager | Predicting near-term changes in the Earth system: A large ensemble of initialized decadal prediction simulations using the Community Earth System Model[END_REF], and the High-Resolution Global Environmental Model (HiGEM [START_REF] Robson | Decadal prediction of the North Atlantic subpolar gyre in the HiGEM high-resolution climate model[END_REF][START_REF] Shaffrey | The new U.K. High-Resolution Global Environment Model: Model description and basic evaluation[END_REF]. Table 1 shows information on the size of the ensembles for the experiments and the resolutions of the models. The models are not necessarily the same as those used in CMIP6 but have been updated since CMIP5 and the external forcings used are those of CMIP5 [START_REF] Taylor | An overview of CMIP5 and the experiment design[END_REF]. The HiGEM model did not participate in CMIP5.

The experiments follow the protocol of Component C of the Decadal Climate Prediction Project (DCPP-C) in CMIP6 [START_REF] Boer | The Decadal Climate Prediction Project (DCPP) contribution to CMIP6[END_REF], they consist of pairs of ensembles, one with all external forcings and one with all external forcings except volcanic aerosols, which are instead set to a background level. The difference between the two (volcano experiment minus no volcano experiment) shows the impacts of the volcanic aerosols. As we only look at the difference between the two states, this method removes any common drifts or offsets from the observed climate. The experiments are initialized with observations with start dates between November and January depending on the model and run for at least 5 years. The volcanic eruption occurs in the first or second year of the experiment (some models do not have hindcasts for every year).

Three volcanoes are included in this study; they are the largest tropical eruptions included in the hindcast period of decadal climate prediction systems: Agung (started erupting March 1963), El Chichón (April 1982), and Pinatubo (June 1991). [START_REF] Stevenson | Role of eruption season in reconciling model and proxy responses to tropical volcanism[END_REF] found that the season the volcanoes erupt in is important for the response. All the volcanoes in this study erupted at a similar time of year, spring to early summer. We average over all three volcanoes taking June of the eruption year to be Month 1. Combining all volcanoes and members results in an ensemble of 135 members with volcanic aerosols and 128 members without volcanic aerosols. This is necessary as mentioned in section 1, a very large ensemble may be needed to capture the model response, which is very likely smaller than in the real world [START_REF] Siegert | A Bayesian framework for verification and recalibration of ensemble forecasts: How uncertain is NAO predictability[END_REF][START_REF] Smith | Robust skill of decadal climate predictions[END_REF].

In four of five models, volcanic aerosols are represented by monthly mean zonal-mean aerosol optical depths (AOD) that are prescribed with a vertical profile in the stratosphere used in the radiation code of the model. EC-Earth, HadGEM3, and HiGEM have 24, 4, and 4 latitude bands, respectively, in which the AOD is constant. The AODs are based on [START_REF] Sato | Stratospheric aerosol optical depths, 1850-1990[END_REF]. CESM has three latitudinal bands, and the AODs are based on [START_REF] Ammann | A monthly and latitudinally varying volcanic forcing dataset in simulations of 20th century climate[END_REF]. MPI-ESM uses zonal-mean monthly mean values of single scattering albedo, aerosol extinction, and asymmetry factor as a function of wavelength and pressure. The data are an extension of [START_REF] Stenchikov | Radiative forcing from the 1991 Mount Pinatubo volcanic eruption[END_REF]. As the implementations of the volcanic forcing are so similar, we do not anticipate the differences to be important, especially against a background of large internal climate variability.

All models used here have slowly varying ozone concentrations, as such, any impacts arising from ozone depletion by the volcanic aerosols will not be represented in these experiments [START_REF] Stenchikov | Arctic Oscillation response to the 1991 Mount Pinatubo eruption: Effects of volcanic aerosols and ozone depletion[END_REF].

Surface temperature observations used in this paper are an average of three observational data sets: The Hadley Centre and Climate Research Unit Surface Temperature Version 4 (HadCRUT4, [START_REF] Morice | Quantifying uncertainties in global and regional temperature change using an ensemble of observational estimates: The HadCRUT4 data set[END_REF], the National Aeronautics and Space Administration-Goddard Institute for Space Studies surface temperature (NASA-GISS, [START_REF] Hansen | Global surface temperature change[END_REF], and the National Climatic Data Center surface temperature (NCDC, [START_REF] Karl | Possible artifacts of data biases in the recent global surface warming hiatus[END_REF]. Air temperature, zonal wind, and mean sea-level pressure observations are an average of three reanalyses: the ECMWF Re-Analysis (ERA40, [START_REF] Uppala | The ERA-40 re-analysis[END_REF], the reanalysis of the National Center for Environmental Protection and National Center for Atmospheric Research (NCEP/NCAR, [START_REF] Kalnay | The NCEP/NCAR 40-year reanalysis project[END_REF], and the Japan Re-Analysis (JRA-55, [START_REF] Harada | The JRA-55 reanalysis: Representation of atmospheric circulation and climate variability[END_REF]. Precipitation observations are taken from the Global Precipitation Climatology Centre (GPCC, [START_REF] Schneider | GPCC's new land surface precipitation climatology based on quality-controlled in situ data and its role in quantifying the global water cycle[END_REF].

Obviously, for observations, there is no equivalent to the experiment without volcanic aerosols, the nonvolcanic reference state has to be approximated in a different way. We try to construct the observed anomalies as analogously as possible to those from the models by similarly making two ensembles, with and without volcanic aerosols. One complicating matter is that the three eruptions studied here all happened during a developing El Niño. Therefore, we decided to only pick initial years for the nonvolcanic reference state that similarly have developing El Niños using the ERSST v5-based data set available from the NCEP website (https://origin.cpc.ncep.noaa.gov/products/analysis_monitoring/ensostuff/ONI_v5.php, last accessed 11 February 2020). We make a nonvolcanic ensemble with six members for the observations by using the four consecutive years following each May of 1972May of , 1976May of , 1977May of , 1979May of , 1986May of , and 1997. The ensemble constructed from these years is subtracted from the three-member ensemble made from the volcanic years starting in 1963, 1982, and 1991 to make the observed anomalies.

For both model and observations, differences are significance tested. Points on maps are considered significant at p ≤ 0.05, and time series have 90% confidence intervals. These are found from the distribution of differences created by bootstrap with replacement of ensemble members from 1,000 tries.

Results

First, we consider the effect of volcanic aerosols on global mean surface temperature to measure the overall climate forcing. Figure 1 shows the multimodel ensemble mean temperature anomaly (blue line) averaged over all three volcanic eruptions (Agung, El Chichón, and Pinatubo). As expected, global mean temperature decreases steadily after the eruptions to reach -0.2 K after about a year. It stays at this level for at least another 18 months.

Figure 1 shows that the response in the observations is much less certain than that of the models. While the global cooling response is barely significant in observations, the response in the models is clearly less than zero apart from the first few months. Estimates made from observations have a large component of noise as many other processes, and forcings were changing at the time of the eruptions. Despite the differences in method and noise, the observations show a similar development and peak cooling to the models. Although the initial cooling appears stronger in the models, there is no significant difference between the monthly means at any time. This implies, as discussed in [START_REF] Lehner | The importance of ENSO phase during volcanic eruptions for detection and attribution[END_REF], that models do not cool too strongly in response to the volcanoes studied here. The result is dependent on the nonvolcanic reference state used to make the observational anomalies also including a developing El Niño, as detailed in section 2. Including or excluding the large El Niño of 1997/1998 in the reference state does not change this result (not shown). Despite being initialized, ENSO has low predictability at a years lead time, so the models will not have captured the observed ENSO phase.

Winds and Atmospheric Temperatures

Figure 2 shows the zonal-mean temperature in the atmosphere throughout the first 2 years after the eruption for the June-July-August (JJA, boreal summer) and December-January-February (DJF, boreal winter)

Figure 2. Zonal-mean cross sections in the atmosphere of temperature anomalies (shaded) and zonal wind anomalies (contours, m/s) for the first and second December-January-February (DJF) and June-July-August (JJA) seasons after the eruptions for the multimodel ensemble mean (a, c, e, g) and an average of reanalyses (b, d, f, h). Values not significant (p > 0.05) for temperature are hatched.

seasons. In general, the most of the lower stratosphere warms and the troposphere cools (particularly in the tropics). The warming in the models and the observations is similar in magnitude, but the tropospheric cooling is patchier and not significant in the observations.

Figure 3a shows a time series of the stratospheric warming in the models (red line) defined as the average temperature at 50 hPa between 30 • S and 30 • N. The peak warming in the stratosphere happens about 7 to 9 months after the eruptions (around the first DJF, DJF1). It then reduces over the next year as the amount of volcanic aerosol declines in the stratosphere. The red line in Figure 3b shows the equivalent time series for the observations. Note that the vertical axis has changed between panels (a) and (b). The timing and magnitude of the stratospheric warming is close between models and observations.

The time evolution in the models of the polar vortex anomalies defined as the average zonal velocity over 50-100 hPa and 55 • N-75 • N is shown by the green line in Figure 3a. The anomalies increase in magnitude in the boreal autumn and peak in late boreal winter in the first year after the eruption. The polar vortex does not exist in boreal summer, but the second DJF (DJF2), the tropical stratospheric heating is a third of what it was in the first winter, and there are no significant anomalies in the polar vortex. The green line in Figure 3b shows that polar vortex on average strengthens much more in the observations than in the models.

The zonal-mean zonal winds (black contour lines) show the strongest response in DJF1, Figures 2c and2d, when the stratospheric warming is the strongest. Even though the stratospheric warming is similar between models and the observations, the response of the wind is stronger in the observations (the polar vortex anomalies for DJF1 in the models and observations are both significantly different from zero, and they are significantly different from each other, defined as that the bootstrapped 90% confidence intervals around the two means do not overlap). The observations are about three times stronger than the models (the range is between 2.5 and 6 times stronger when taking the uncertainty of the model mean into account).

It appears from Figures 2c and2d that the northern edge of the midlatitude tropospheric jet is weakened in DJF1. The latitude of the jet changes with longitude, so the zonal-mean value is misleading. In the Atlantic, the jet is split into two: The midlatitude jet is shifted North, while the subtropical jet is strengthened. In the Pacific, the jet is shifted South (see in supporting information Figure S1).

10.1029/2019JD031739 Also notable in Figures 2c,2d, and 2h are the cold temperature anomalies at high latitudes in boreal winter. In DJF1, like the zonal winds, the cooling anomaly is much weaker in models compared to observations. The time series of the cooling in the Northern Hemisphere defined as the average temperature over 50-200 hPa and 70 • N-90 • N is shown in Figure 3 as the blue line. For both models and observations, it is anticorrelated with the polar vortex. This implies that the upper-troposphere cold anomalies over the Arctic are probably caused by the strengthened polar vortex through cooling as air rises [START_REF] Limpasuvan | Stratosphere-troposphere evolution during polar vortex intensification[END_REF].

In DJF2, Figures 2g and2h, models show a weak decrease in the polar vortex, while in observations, the response is similar to DJF1, although it is weaker (the polar vortex anomaly for observations is significantly different from zero and significantly different from the models, the polar vortex response for the models is not significantly different from zero). Associated with this, the upper-troposphere Arctic temperatures are anomalously cold in the observations, but there is no significant anomaly in the models.

Surface Responses

The Northern Hemisphere circulation changes in the mean sea-level pressure in the two winters following the eruptions are shown in Figure 4. In DJF1, panels (a) and (b), there is a low over the Arctic and a high over Europe for both models and observations. [START_REF] Driscoll | Coupled Model Intercomparison Project 5 (CMIP5) simulations of climate following volcanic eruptions[END_REF] found that the response in CMIP5 models was about 10 times weaker than in observations. We find a similar result here, the model anomalies are about seven times weaker than the observational anomalies in the Atlantic region. We have used a different range of contours for models and observations in Figure 4 so that patterns can be compared even if the magnitudes are different. As found above, in DJF2, observations have a similar response to DJF1, compare panels (b) and (d). The response in the models, shown in panel (c), is different and is dominated by a low pressure in the Northeast Atlantic. This circulation change in DJF1 has an impact on Northern Hemisphere winter surface temperatures, shown in Figures 5a and5b. While in DJF1, the tropics have cooled (especially over land), both the models and the observations show warm anomalies over Scandinavia and northern Russia (though this is not significant for the models). Again, model anomalies are weaker than in the observations for the surface temperature, presumably due the weaker circulation response. In DJF1 and DJF2, the cooling appears more uniform in the models than in the observations, though this may be because we are comparing an ensemble mean with an average of three volcanoes, which would have more noise. As discussed above, the global mean surface temperature response is not significantly different between the observations and the models.

Tropical precipitation, shown in Figure 6, has anomalies over equatorial Africa already in JJA1 (panels a and b), though they are not significant. This is the largest equatorial land mass, and as it cools faster than the surrounding oceans, this might cause the reduction in precipitation [START_REF] Khodri | Tropical explosive volcanic eruptions can trigger El Niño by cooling tropical Africa[END_REF]. The main feature of the precipitation, for which there is good agreement between models and observations, especially after JJA1, is that precipitation is reduced over the maritime continent. In DJF1 (panel c), the models show positive precipitation anomalies in the equatorial West Pacific. This is maintained for the next year as the precipitation anomalies spread East (panels e and g). There are dry anomalies in Australia in DJF1 (panels c and d) and in Southeast Asia in JJA2 (panels e and f) that suggest the monsoons there are reduced as was found by several studies, for example, [START_REF] Iles | The effect of volcanic eruptions on global precipitation[END_REF], [START_REF] Iles | The global precipitation response to volcanic eruptions in the CMIP5 models[END_REF], [START_REF] Liu | Global monsoon precipitation responses to large volcanic eruptions[END_REF][START_REF] Paik | Assessing the impact of volcanic eruptions on climate extremes using CMIP5 models[END_REF][START_REF] Paik | Assessing the impact of volcanic eruptions on climate extremes using CMIP5 models[END_REF].

Impact on ENSO

The precipitation in Figure 6 shows, at least for the models, indications of the growth of an El Niño-like state from JJA1 to DJF2. Precipitation reduces over the Maritime Continent and increases to the East in the first year after the eruption. This increase in precipitation then spreads eastwards over the following year. For DJF2, the temperatures of Figure 5c in the central equatorial Pacific show no cooling in contrast to most of the rest of the tropics.

To see this more clearly, Figure 7 shows Hovmöller plots of monthly mean model precipitation and temperature anomalies in the equatorial Indo-Pacific. After about 8 months, around January of the first year after the eruptions, a clear dipole emerges in the precipitation, indicating that it has shifted East. In the following year, it spreads further East. As surface temperatures are dominated by the tropical cooling, in Figure 7b, we 10.1029/2019JD031739 show the surface temperature relative to the whole of the tropics (20 • S-20 • N). This shows eastward propagation of anomalies from about a year after the eruption, peaking in the second October after the eruption, similar to an El Niño.

Figure 8 shows relative sea surface temperature (RSST), the temperature difference between the NINO3.4 region and the tropical average (20 • S-20 • N) as used by [START_REF] Khodri | Tropical explosive volcanic eruptions can trigger El Niño by cooling tropical Africa[END_REF]. For the models, RSST is positive for the first 2 years after the eruptions and peaks in the second October. The observed average is also positive or stays close to zero. The uncertainties around the observed values are large, and the confidence interval never excludes zero but also includes the model ensemble.

Sensitivity studies by [START_REF] Khodri | Tropical explosive volcanic eruptions can trigger El Niño by cooling tropical Africa[END_REF] in the atmospheric component of the Institute Pierre Simon Laplace-Coupled Model (IPSL-CM5B, not part of this study) showed that the chain of events leading to an El Niño-like response to volcanic aerosols in the second year started with reduced African precipitation. Figure 6a shows precipitation anomalies of the right sign over Africa in JJA1 (for models only), and Figure 7a shows propagation eastward in the Indian Ocean following this (not significant). So while we cannot exclude Africa being an important initiator of the mechanism, the strongest signals appear over the Maritime Continent about 6 months after the eruptions.

The results here are consistent with [START_REF] Ohba | Impact of strong tropical volcanic eruptions on ENSO simulated in a coupled GCM[END_REF] and [START_REF] Predybaylo | Impacts of a Pinatubo-size volcanic eruption on ENSO[END_REF] who find that the Maritime Continent is important for the ENSO response. It is the relatively stronger cooling here compared to the tropical West Pacific that likely causes the precipitation to move to the tropical West Pacific. The sea 10.1029/2019JD031739 surface temperature (SST) in the central and East Pacific is controlled by upwelling and so cools less than the West Pacific. This weakens the easterly trade winds, which increases the SSTs in the East Pacific further through the Bjerknes feedback, leading to an El Niño response [START_REF] Predybaylo | Impacts of a Pinatubo-size volcanic eruption on ENSO[END_REF].

Despite the lack of a detectable El Niño in observations, the precipitation fields in Figure 6h for observations in DJF2 show many known teleconnections to El Niño: low precipitation in Southern Africa, Australia, and the Maritime Continent as well as increased precipitation in Southern United States. [START_REF] Stevenson | El Niño like" hydroclimate responses to last millennium volcanic eruptions[END_REF] Figure 8. Surface temperature relative to the tropics in the NINO3.4 region for the multimodel ensemble (blue) and observations (black). The shaded areas around the time series represent the 90% confidence intervals found from a bootstrap with replacement. found that circulation changes due to volcanic eruptions can lead to precipitation anomalies that are similar to those associated with El Niño even when there is no ENSO event in surface temperatures.

Impact on the NAO

We now look closer at the climate response over the Atlantic sector. Figure 9 shows the model NAO index and observed NAO. The NAO is defined as the difference, south minus north, in mean sea-level pressure between two large boxes bounded by the longitudes 90 • W-60 • E. The northern box is bounded by 55 • N-90 • N and the southern box by 20 • N-55 • N as in [START_REF] Stephenson | North Atlantic Oscillation response to transient greenhouse gas forcing and the impact on European winter climate: A CMIP2 multi-model assessment[END_REF]. The boxes are shown in Figure 4a and are chosen to be large so that the exact location of the NAO in the different models is not important.

The NAO indices in the model and observations are both positive in the first winter after the eruption, but as expected from Figures 4a and4b, the amplitude is about seven times stronger in the observations than in the models (note that the 90% confidence intervals do not overlap between models and observations for the first January). In the second winter, the response is not significant as the confidence intervals for the observed NAO and modeled NAO both include zero.

Comparing the NAO (Figure 9) with the indices in Figure 3 supports that the NAO response in DJF1 is related to the stratospheric tropical warming via the polar vortex. The warming increases the stratospheric temperature gradient between 40 • N and the pole. This strengthens the polar vortex [START_REF] Stenchikov | Arctic Oscillation response to the 1991 Mount Pinatubo eruption: Effects of volcanic aerosols and ozone depletion[END_REF][START_REF] Shindell | Dynamic winter climate response to large tropical volcanic eruptions since 1600[END_REF][START_REF] Driscoll | Coupled Model Intercomparison Project 5 (CMIP5) simulations of climate following volcanic eruptions[END_REF][START_REF] Swingedouw | Impact of explosive volcanic eruptions on the main climate variability modes[END_REF]. In addition, increasing extratropical westerly wind anomalies (Figures 2c and2d) increase propagation of planetary waves into the tropics and hence reduce the wave flux into the vortex, further stabilizing and strengthening it (Bittner et al., 2016), though this has not been analyzed here. The lag between the tropical stratospheric warming and the strengthening of the polar vortex at the beginning of the time series in Figure 3 is likely because the polar vortex only exists in winter.

The tropical stratospheric warming has been questioned as a mechanism [START_REF] Swingedouw | Impact of explosive volcanic eruptions on the main climate variability modes[END_REF] as in the stratosphere, the tropics are colder than the polar region. This would mean that warming the tropical stratosphere would decrease the meridional temperature gradient, not increase it. However, in the models used here, we find in DJF1 that there is a climatological temperature maximum (found by averaging all DJF1 of the experiments without volcanic aerosols) at 40 • N at the 50-and 100-hPa levels, so the warming does indeed increase the temperature gradient. This climatological maximum exists in the observations as well, so this is a realistic response.

Impact on the North Atlantic Ocean

To explore the potential longer-term climate response (though the experiments only continue for 4 years after the volcanic eruption), we look at the North Atlantic Ocean, a region where volcanoes can induce changes at decadal timescales [START_REF] Ding | Ocean response to volcanic eruptions in Coupled Model Intercomparison Project 5 simulations[END_REF][START_REF] Mignot | Volcanic impact on the Atlantic Ocean over the last millennium[END_REF][START_REF] Stenchikov | Volcanic signals in oceans[END_REF][START_REF] Swingedouw | Impact of explosive volcanic eruptions on the main climate variability modes[END_REF]. Figure 10a shows an anomalously deep mixed layer (which can be interpreted as a proxy for deep convection) during the three boreal winters following the eruptions, with a maximum in the Labrador Sea. The time series averaged over the blue box can be seen in Figure 10b. It shows that the ensemble with volcanic aerosols included has a deeper mixed layer than the ensemble without volcanic aerosols in late boreal winter and spring for the first 2 years after the eruptions. A strengthened NAO, as seen in Figure 9 for the first winter and early spring, would favor downwelling from Ekman pumping south of the strengthened winds (where there is negative wind stress curl), leading to the increased mixed layer depth (MLD) seen in Figure 10 for the first winter.

Figure 10c shows the AMOC streamfunction anomaly for Years 2-4 after the volcanic eruption. A barotropic response can be seen, strongest between 30 • N and 40 • N, which strengthens the circulation. Although the anomaly is not significant over this period, this is in line with previous literature identifying similar AMOC intensifications in different models, which are caused by either a rapid dynamical adjustment related to a dominant negative wind stress curl anomaly over the subpolar North Atlantic forced by the eruptions [START_REF] Mignot | Volcanic impact on the Atlantic Ocean over the last millennium[END_REF][START_REF] Zanchettin | Bi-decadal variability excited in the coupled ocean-atmosphere system by strong tropical volcanic eruptions[END_REF] or enhanced deep convection activity in response to the surface cooling, salinification (due to reduced precipitation), and subsequent weakening of density stratification in the Labrador Sea [START_REF] Iwi | Mechanisms linking volcanic aerosols to the Atlantic Meridional Overturning Circulation[END_REF][START_REF] Ortega | Variability of the Atlantic meridional overturning circulation in the last millennium and two IPCC scenarios[END_REF][START_REF] Stenchikov | Volcanic signals in oceans[END_REF]. There is a significant response in the streamfunction over the equator, which is caused by anomalous easterly winds. NAO signal weakens after the first winter, the response to surface cooling and salinification is likely to become more important in driving the AMOC anomalies thereafter [START_REF] Stenchikov | Volcanic signals in oceans[END_REF].

Conclusions

In this multimodel study with initialized climate prediction systems, we have used an experimental design created specifically to clarify the impact of volcanic aerosols on near-term climate and reduce model biases [START_REF] Boer | The Decadal Climate Prediction Project (DCPP) contribution to CMIP6[END_REF]: Anomalies are the difference between an ensemble that includes volcanic aerosols and one without volcanic aerosols. Both ensembles include all other external forcings relevant to the time of the eruptions. We combine three different volcanic eruptions, Agung, El Chichón, and Pinatubo, to get a total ensemble size of 135 members with volcanic aerosols and 128 without.

Global temperature cools, peaking at about -0.2 K a year after the eruptions in both models and observations (Figure 1), though uncertainties are large in the observations. We find that the models simulate realistic levels of warming in the stratosphere (Figure 2), but zonal-mean winds in the first winter are approximately three times weaker than in observations (cf. Figures 2c and2d). At the surface, both mean sea-level pressure (Figures 4a and4b) and temperature (Figures 5a and5b) show similar patterns to observations in the first winter, but the modeled amplitude is about seven times weaker. We are not claiming that models cannot display zonal winds, temperatures, and sea-level pressures of the observed magnitude, but since the model response is weak in the ensemble mean (which is a measure of the signal, after the noise has been averaged out), then either the response is spurious in the observations or the response is too weak in the models.

Comparison with observations is difficult as they are noisy, and we are only considering a few cases, but this makes the similarity in patterns between models and observations all the more surprising and implies that the observed signal is stronger than previously thought. The case for the response being too weak in the models is strengthened by the fact that many recent studies have found that patterns in supposedly noisy observations are reproduced in models, but with reduced magnitude, especially around the North Atlantic [START_REF] Baker | An intercomparison of skill and overconfidence/underconfidence of the wintertime North Atlantic Oscillation in multimodel seasonal forecasts[END_REF][START_REF] Dunstone | Skilful predictions of the winter North Atlantic Oscillation one year ahead[END_REF][START_REF] Scaife | A signal-to-noise paradox in climate science[END_REF][START_REF] Siegert | A Bayesian framework for verification and recalibration of ensemble forecasts: How uncertain is NAO predictability[END_REF][START_REF] Smith | Robust skill of decadal climate predictions[END_REF]. In the second winter, observations show a persistence of the same pattern as the first winter, but it is weaker and less significant. This persistence is not simulated by the models.

An El Niño-like state, which develops in the models about 1 year after the eruptions, grows to a maximum around the second winter (Figure 7). This is initiated by an eastward shift in precipitation from the Maritime Continent toward the western tropical Pacific that starts within a few months of the eruptions (Figures 6a and6b). This could be related to the Maritime Continent cooling faster under the influence of volcanic aerosols than the surrounding oceans. This robust impact on ENSO in the models is not visible in the observations, possibly due to the large amount of noise present (Figure 8).

Models and observations both have a positive NAO-like pattern in the first winter after the eruptions (Figures 4a and4b), in agreement with some previous studies [START_REF] Shindell | Dynamic winter climate response to large tropical volcanic eruptions since 1600[END_REF][START_REF] Stenchikov | Arctic Oscillation response to the 1991 Mount Pinatubo eruption: Effects of volcanic aerosols and ozone depletion[END_REF][START_REF] Zambri | Winter warming and summer monsoon reduction after volcanic eruptions in Coupled Model Intercomparison Project 5 (CMIP5) simulations[END_REF]. The positive NAO is likely caused by a strengthened polar vortex driven by the warming in the lower tropical stratosphere (Figure 3). The amplitude of the NAO anomaly is about seven times weaker in the models than in the observations in the first winter (Figure 9). The signal-to-noise ratio appears to be too small in the models, consistent with other studies [START_REF] Baker | An intercomparison of skill and overconfidence/underconfidence of the wintertime North Atlantic Oscillation in multimodel seasonal forecasts[END_REF][START_REF] Dunstone | Skilful predictions of the winter North Atlantic Oscillation one year ahead[END_REF][START_REF] Eade | Do seasonal-to-decadal climate predictions underestimate the predictability of the real world[END_REF][START_REF] Scaife | A signal-to-noise paradox in climate science[END_REF]. In the second winter, neither the observations nor the models have a significant NAO signal, though the observed mean sea-level pattern is NAO-like (Figure 4d).

Finally, we find dynamical changes in the Atlantic Ocean (Figure 10). Even for the short length of these experiments (4 years), there is already a strengthening of the AMOC. This is probably mediated by changes in wind stress curl and increased MLDs. The volcanic response therefore likely extends beyond the length of these experiments, potentially leading to a positive phase of the Atlantic Multidecadal Variability (AMV, [START_REF] Knight | A signature of persistent natural thermohaline circulation cycles in observed climate[END_REF][START_REF] Ottera | External forcing as a metronome for Atlantic multidecadal variability[END_REF] and increases to the ocean heat content in the subpolar North Atlantic [START_REF] Hermanson | Forecast cooling of the Atlantic subpolar gyre and associated impacts[END_REF][START_REF] Robson | Causes of the rapid warming of the North Atlantic Ocean in the mid-1990s[END_REF], with important impacts over land [START_REF] Hermanson | Forecast cooling of the Atlantic subpolar gyre and associated impacts[END_REF][START_REF] Robson | Predictable climate impacts of the decadal changes in the ocean in the 1990s[END_REF].

Figure 1 .

 1 Figure 1.Global temperature changes due to volcanic aerosols for the multimodel ensemble (blue) and observations (black). The numbers on the month labels on the time axis refer to the number of years since the eruption. The shaded areas around the time series represent the 90% confidence intervals found from a bootstrap with replacement.

Figure 3 .

 3 Figure 3. Zonal wind and temperature anomalies for the multimodel ensemble (a) and observations (b) with 90% confidence intervals. Polar vortex anomalies (PolarVortex, green) are the average zonal velocity over 50-100 hPa and 55 • N-75 • N, high-latitude upper-troposphere temperature (HiLatitudeT, blue) is the average temperature over 50-200 hPa and 70 • N-90 • N, and stratospheric temperature (StratosphericT, red) is the average temperature at 50 hPa between 30 • S and 30 • N. Note that the vertical axis scale is different between (a) and (b).

Figure 4 .

 4 Figure 4. Mean sea-level pressure anomalies for the first (a, b) and second DJF (c, d) after the eruptions for the multimodel ensemble mean (a, c) and an average of reanalyses (b, d). Note that the color scales are different between models and observations. Significant values are stippled.

Figure 5 .

 5 Figure 5. Surface temperature anomalies for the first (a, b) and second DJF (c, d) after the eruptions for the multimodel ensemble mean (a, c) and an average of reanalyses (b, d). Significant values are stippled.

Figure 6 .

 6 Figure 6. Precipitation anomalies (mm/day) for the first and second DJF and JJA after the eruptions for the multimodel ensemble mean (a, c, e, g) and observations (b, d, f, h). Significant values are stippled. Note that observations do not exist over the ocean (gray).

Figure 7 .

 7 Figure 7. Indo-Pacific equatorial (5 • S-5 • N) multimodel ensemble mean Hovmöller plots for (a) precipitation (mm/day) and (b) surface temperature relative to the tropics (20 • S-20 • N, K). Stippling indicates significant values.

Figure 9 .

 9 Figure 9. NAO index for the multimodel ensemble (blue) and observations (black). The shaded areas around the time series represent the 90% confidence intervals found from a bootstrap with replacement.

Figure 10 .

 10 Figure 10. North Atlantic dynamic impacts. (a) Multimodel ensemble mean mixed layer depth (MLD) anomaly averaged over the first three winters (DJF) following the volcanic eruptions and (b) averaged monthly MLD anomaly over the Labrador Sea (blue box in panel a). (c) Multimodel ensemble mean Atlantic Meridional Overturning Circulation (AMOC) streamfunction anomaly averaged over Years 2-4 after the eruptions and (d) averaged AMOC streamfunction anomaly (green box in panel c). A 12-month running mean has been applied to this time series to focus on the interannual variability.

A

  smoothed time series of the AMOC from the green box on Figure 10c is shown in Figure 10d. The response of the AMOC to the deepened MLD typically has a lag of 2-3 years between a peak in convection and the AMOC (Eden & Willebrand, 2001; Ortega et al., 2017), which is consistent with Figure 10d. As the positive 10.1029/2019JD031739

Table 1

 1

  Table of the Participating Models

	Model	Ensemble size	Resolution
	EC-Earth V2.3	5	T159 L62, 1 • L42
	MPI-ESM	10 (3)	T63 L47, 1.5 • L40
	HadGEM3-GC2	10	N216 L85, 0.25 • L75
	CESM 1.1	10	1 • L30, 1 • L60
	HiGEM	10	0.83 • × 1.25 • L38, 1/3 • L40
	Note. The ensemble size refers to the number of ensemble members per
	experiment (counting the with and without volcanic aerosols experi-
	ments as separate). MPI-ESM has 10 members in its ensemble for most
	experiments, except for the Agung experiment without volcanic aerosols,
	which has three. Resolution is described with the atmospheric resolu-
	tion first and then oceanic resolution, and the number following letter L
	indicates the number of vertical levels.

Acknowledgments

We thank the three reviewers for their help in improving this manuscript. The data are available from DOI (10.5281/zenodo.2613699). L. H. was

Adams, J. B., Mann, M. E., & Ammann, C. M. (2003). Proxy evidence for an El Niño-like response to volcanic forcing. Nature, 426, 274-278. https://doi.org/10.1038/nature02101