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 (Ding et al., 2014;Robock, 2000;Timmreck, 2012;Timmreck et al., 2016;Zhong et al., 2011).Volcanic impacts are needed to fully explain the climate variability of the twentieth century (Smith et al., 2016).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 (Hansen et al., 1997;Lacis et al., 1992;Swingedouw et al., 2017).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 (Stenchikov et al., 2002).The result is an increased chance of a positive North Atlantic Oscillation (NAO) after the volcanic eruption, as discussed by many authors including Robock and Mao (1992), Robock (2000), Stenchikov et al. (2002), Shindell et al. (2004), and Stenchikov et al. (2006).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 (Dunstone et al., 2016).Ding et al. (2014) 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 (Khodri et al., 2017;Predybaylo et al., 2017;Stevenson et al., 2016).There 10.1029/2019JD031739 is also compelling evidence of this from paleoreconstructions (Adams et al., 2003;Wilson et al., 2010).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 (Swingedouw et al., 2017).
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 (Iwi et al., 2012;Ottera et al., 2010;Stenchikov et al., 2009;Swingedouw et al., 2017).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 (Ding et al., 2014), though changes in windstress curl are also likely to play a role (Mignot et al., 2011).In some models, background conditions can alter the response in the North Atlantic (Swingedouw et al., 2015;Zanchettin et al., 2013).
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., Driscoll et al., 2012).However, experiment design is important, the Pacific response may depend on capturing the phase of ENSO at the time of the eruption (Lehner et al., 2016).Commonly used seasonal averages may not capture the response (Barnes et al., 2016;Predybaylo et al., 2017).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 (Baker et al., 2018;Dunstone et al., 2016;Scaife & Smith, 2018).
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;Timmreck et al., 2016).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 (Illing et al., 2018) 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;Driscoll et al., 2012;Iles & Hegerl, 2014;Khodri et al., 2017;Maher et al., 2015;Zambri & Robock, 2016).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 (Clark et al., 2017).

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 (Hazeleger et al., 2012;Ménégoz et al., 2018), the Max Planck Institute Earth System Model (MPI-ESM, Giorgetta et al., 2013;Pohlmann et al., 2013), the Hadley Centre Global Environmental Model Global Coupled Configuration 2 (HadGEM3-GC2, Dunstone et al., 2016, Williams et al., 2015), the Community Earth System Model (CESM) 1.1 (Hurrell et al., 2013;Yeager et al., 2018), and the High-Resolution Global Environmental Model (HiGEM Robson et al., 2018;Shaffrey et al., 2009).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 (Taylor et al., 2012).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 (Boer et al., 2016), 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 10.1029/2019JD031739 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).Stevenson et al. (2017) 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 (Siegert et al., 2016;Smith et al., 2019).
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 Sato et al. (1993).CESM has three latitudinal bands, and the AODs are based on Ammann et al. (2003).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 Stenchikov et al. (1998).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 (Stenchikov et al., 2002).
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, Morice et al., 2012), the National Aeronautics and Space Administration-Goddard Institute for Space Studies surface temperature (NASA-GISS, Hansen et al., 2010), and the National Climatic Data Center surface temperature (NCDC, Karl et al., 2015).Air temperature, zonal wind, and mean sea-level pressure observations are an average of three reanalyses: the ECMWF Re-Analysis (ERA40, Uppala et al., 2005), the reanalysis of the National Center for Environmental Protection and National Center for Atmospheric Research (NCEP/NCAR, Kalnay et al., 1996), and the Japan Re-Analysis (JRA-55, Harada et al., 2016).Precipitation observations are taken from the Global Precipitation Climatology Centre (GPCC, Schneider et al., 2014).
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 Lehner et al. (2016), 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 and 2d, 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 and 2d 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).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 (Limpasuvan et al., 2005).
In DJF2, Figures 2g and 2h, 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.Driscoll et al. (2012) 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 and 5b.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 (Khodri et al., 2017).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, Iles et al. (2013), Iles and Hegerl (2014), Liu et al. (2016), andPaik andMin (2018).

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.4region and the tropical average (20 • S-20 • N) as used by Khodri et al. (2017).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 Khodri et al. (2017) 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 Ohba et al. (2013) and Predybaylo et al. (2017) 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 (Predybaylo et al., 2017).
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.Stevenson et al. (2016) Figure 8. Surface temperature relative to the tropics in the NINO3.4region 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 Stephenson et al. (2006).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 and 4b, 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 (Stenchikov et al., 2002;Shindell et al., 2004;Driscoll et al., 2012;Swingedouw et al., 2017).In addition, increasing extratropical westerly wind anomalies (Figures 2c and 2d) 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 (Swingedouw et al., 2017) 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 10.1029/2019JD031739 changes at decadal timescales (Ding et al., 2014;Mignot et al., 2011;Stenchikov et al., 2009;Swingedouw et al., 2017).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 (Mignot et al., 2011;Zanchettin et al., 2012) 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 (Iwi et al., 2012;Ortega et al., 2012;Stenchikov et al., 2009).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 (Stenchikov et al., 2009).

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 (Boer et al., 2016): 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 and 2d).At the surface, both mean sea-level pressure (Figures 4a and 4b) and temperature (Figures 5a and 5b) 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 (Baker et al., 2018;Dunstone et al., 2016;Scaife & Smith, 2018;Siegert et al., 2016;Smith et al., 2019).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 and 6b).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 and 4b), in agreement with some previous studies (Shindell et al., 2004;Stenchikov et al., 2002;Zambri & Robock, 2016).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 (Baker et al., 2018;Dunstone et al., 2016;Eade et al., 2014;Scaife & Smith, 2018).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, Knight et al., 2005;Ottera et al., 2010) and increases to the ocean heat content in the subpolar North Atlantic (Hermanson et al., 2014;Robson et al., 2012), with important impacts over land (Hermanson et al., 2014;Robson et al., 2013).

Figure 1 .
Figure1.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 .
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 .
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 .
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 .
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 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 .
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
Table of the Participating Models