The impact of ENSO and NAO initial conditions and anomalies on the modeled response to Pinatubo-sized volcanic forcing

Strong, strato-volcanic eruptions are a substantial, intermittent source of natural climate variability. Initial atmospheric and oceanic conditions, such as El Niño Southern Oscillation (ENSO) and the North Atlantic Oscillation (NAO), also naturally impact climate on interannual timescales. We examine how initial conditions of ENSO and NAO contribute to the evolution of climate in the period following a Pinatubo-type eruption using a large (81-member) ensemble of model simulations in GISS model E2.1-G. Simulations are initialized from sampled conditions of ENSO and NAO using the protocol of the coordinated CMIP6 Volcanic Model Intercomparison Project (VolMIP) – where aerosols are forced with respect to time, latitude, and height. We analyze paired anomalous variations (perturbed – control) to understand changes in global and regional climate responses under positive, negative, and neutral ENSO and NAO conditions. In particular, we find that for paired anomalies there is a high probability of strong (∼1.5 ^∘C) warming of northern Eurasia surface air temperature in the first winter after the volcanic eruption for negative NAO ensembles coincident with decreased lower stratospheric temperature at the poles, decreased geopotential height, and strengthening of the stratospheric polar vortex. Climate anomalies (relative to average conditions across the control period), however, show no mean warming and suggest that the strength of this response is impacted by conditions present in the selected period of the control run. Again using paired anomalies, we also observe that under both +ENSO and −ENSO ensembles sea surface temperature decreases in the first post-eruptive boreal winter coinciding with surface cooling from volcanic aerosols. Neutral ENSO ensembles, on the other hand, show variability in their response with no clear trend in post-eruptive warming or cooling. In general, paired anomalies from unperturbed simulations give insight into the evolution of the climate response to volcanic forcing; however, when compared with anomalies from climatological conditions, it is clear that paired anomalies are significantly affected by sampled initial conditions occurring at the time of the volcanic eruption.

1 Introduction

Strong, explosive volcanic eruptions are an intermittent, natural source of climate variability acting on both interannual and decadal scales. Explosive volcanic eruptions eject sulfur dioxide, halogens, ash, and water vapor into the stratosphere, where the particles are converted into sulfate aerosols (LeGrande et al., 2016). The loading of stratospheric aerosols increases aerosol optical depth of the atmosphere (Lacis, 2015), thus imposing a radiative forcing via scattering of shortwave radiation and absorption of longwave radiation in the stratosphere (Zanchettin et al., 2013).

The impact that a strong volcanic eruption has on the climate system depends on many factors including size, ejection height, and location of the eruption. These factors impact the amount, location, and dynamics of how aerosols are loaded in the atmosphere (Timmreck et al., 2010; LeGrande and Anchukaitis, 2015). Mount Pinatubo is an example of one such strong volcanic eruption which occurred in the Philippines in June 1991, ejecting 18 Tg of SO₂ into the atmosphere at a height of 20 km (Stenchikov et al., 1998). The Pinatubo eruption is widely studied as it is one of the largest volcanic eruptions since 1991, thus providing well-constrained satellite observations to run climate model simulations (McCormick et al., 1995; Bluth et al., 1992; Stenchikov et al., 1998).

The timing of an eruption, including the time of season and initial climate conditions, also impact climate's response to a volcanic eruption by actively influencing the mechanisms involved in the post-eruption decadal climate evolution (Zanchettin et al., 2013). Climate variability, such as the El Niño Southern Oscillation System (ENSO) and North Atlantic Oscillation (NAO), continuously cause variations in earth's climate over time (Philander, 1983; Allan et al., 1996; Timmermann et al., 2018; Hurrell et al., 2003; Wanner et al., 2001). Climate modeling studies have examined how the initial state of these climate conditions impact the climate response to volcanic eruption on interannual to decadal scales: while radiative forcing impacts remained the same, different initial climate states caused substantial variability in surface atmospheric and oceanic conditions (Zanchettin et al., 2013; Pausata et al., 2020). Further climate modeling experiments have thus been designed to capture variability that may occur due to different initial states of the climate system at the time of a modeled volcanic eruption (Zanchettin et al., 2016).

In this study, we analyze the modeled response to a 1991 Pinatubo-sized eruption under the coordinated protocol of the Volcanic Model Intercomparison Project (VolMIP) (Eyring et al., 2016) aiming to understand how initial ENSO and NAO climate conditions can impact the climate response. In particular, we investigate how initial ENSO and NAO conditions, which are sampled for as part of the VolMIP protocol (Zanchettin et al., 2016), impact the response in the ENSO region and in the Northern Hemisphere for the first 3 years after the eruption. We also use this ensemble to examine the impact of initial climate conditions and initial climate sampling on the aggregate response of the climate conditions to volcanic forcing.

2 Relevant volcanic climate responses

2.1 Impact of the Pinatubo eruption on climate

The impact of a Pinatubo-sized eruption on the earth's climate system is significant. Previous work has shown that Pinatubo-sized volcanic eruptions decrease radiative flux in the region (40^∘ N–40^∘ S) by around −4.3 W m⁻² at their peak aerosol forcing (Minnis et al., 1993), with radiative effects lasting for about 2 years after the eruption. In comparison, anthropogenic radiative forcing is estimated to have increased global energy budget by 2.3 W m⁻² over the industrial period (Myhre et al., 2013), making volcanic forcing a short-lived but substantial source of natural climate variability. The resulting impacts in the climate system, however, last years after volcanic aerosols have been depleted. The direct impacts of volcanic aerosols include cooling of the earth's surface and warming of the stratosphere (Lacis, 2015). These direct impacts initiate many other changes in the climate system including changes in atmospheric circulation, the hydrological cycle, the cryosphere, and the carbon cycle (Zanchettin et al., 2016). Below we discuss specifically two responses: the response of ENSO and the response of the climate system in the Northern Hemisphere.

2.2 ENSO response

El Niño Southern Oscillation (ENSO) is an important mode of climate variability which oscillates between positive (El Niño), neutral, and negative (La Niña) phases at timescales of about 2–7 years in the equatorial Pacific Ocean (Predybaylo et al., 2017). During positive phases, the equatorial Pacific experiences higher than average sea surface temperature, and during negative phases, the average sea surface temperature is lower. These oceanic changes are associated with changes in both regional climate and global climate connections. Both observational (direct and proxy-based) and model-based studies have been used to examine ENSO responses to large volcanic perturbations. Some proxy-based and several modeling studies suggest that large, tropical volcanic eruptions increase the likelihood of an El Niño-like, or positive, sea surface temperature anomaly following the eruption (Adams et al., 2003; Predybaylo et al., 2017; Khodri et al., 2017).

Predybaylo et al. (2017) and Zambri et al. (2019) additionally studied the robustness of the simulated El Niño anomaly under varying initial conditions at the time of volcanic eruptions. While Predybaylo et al. (2017) found enhanced El Niño-like warming for all Pinatubo simulations except those where eruptions occurred in La Niña years, Zambri et al. (2019) found a consistent warming of tropical sea surface temperature of 0.5–1.0 ^∘C in the Niño 3.4 region in response to the 1783 Laki Eruption in the Whole Atmosphere Community Climate Model (WACCM). Research has also focused on understanding the dynamics of the El Niño anomaly. Suggested mechanisms for this increase in sea surface temperature includes the ocean dynamical thermostat (Clement et al., 1996), which causes advection of warm water through differential cooling, and Walker circulation cells which are caused by post-eruptive cooling over tropical Africa (Khodri et al., 2017). This mechanism was also shown to cause sustained 7-year El Niño anomalies in response to soot aerosols from simulated global nuclear war (Coupe et al., 2021). Despite several studies supporting El Niño like anomalies, other observational and modeling studies suggest that there is no statistically significant El Niño-like response after several large volcanic eruptions (Dee et al., 2020). These studies suggest that anomalies found in observational records and model simulations are not statistically significant and are rather within the range of natural climate variability (Dee et al., 2020).

2.3 Northern Hemisphere winter response

The Northern Hemisphere experiences a unique response during the first winter after large volcanic eruptions. Many observational (Graf et al., 2007; Christiansen, 2008) and modeling (Timmreck, 2012; Stenchikov et al., 2002) studies have noted a strengthening of the polar vortex the first winter after a large volcanic eruption. This increased polar vortex circulation in the lower stratosphere is closely associated with an enhanced phase of the Arctic Oscillation (AO) and North Atlantic Oscillation (NAO) – two modes of natural climate variability that are separately defined but closely related in their associated climate impacts including surface temperature patterns (Cohen and Barlow, 2005). Such increased surface air temperature patterns have commonly been observed after large volcanic eruptions such as the 1991 Pinatubo eruption (Robock and Mao, 1995; Kelly et al., 1996). Thus, the unique signature of increased surface air temperature over Eurasia termed “winter warming” has been analyzed in several volcanic modeling studies.

Modeling studies from previous climate model intercomparison projects substantiate post-eruptive winter warming. For example, Zambri and Robock (2016) analyzed an ensemble of CMIP5 simulations finding that most models produce a winter warming signature over the Northern Hemisphere corresponding with a stronger polar vortex in the lower stratosphere both over the historical 1850–2005 simulation period (Zambri and Robock, 2016) and over the past millennium (Zambri et al., 2017). Analyses from individual models have also previously supported winter warming corresponding with strengthened polar vortex circulation. For example, the NCAR CAM5 AMIP large ensemble showed consistent winter warming in response to both the 1982 El Chinchón and the 1991 Pinatubo eruptions (Coupe and Robock, 2021). This increase in surface temperature is also seen in both observational and global modeling studies (Robock and Mao, 1992; Graft et al., 1993).

Still other studies call the robustness of this modeled result into question. For example, other analyses of CMIP5 models show variation in the prevalence of this response (Timmreck et al., 2016; Driscoll et al., 2012) suggesting that large numbers of ensembles may be required to see a significant strengthening of the polar vortex (Bittner et al., 2016). One proposed cause for inconsistencies in the winter warming response is that the simulated winter warming response in a model is within the range of internal variability (Polvani et al., 2019) and thus is not a robust response to volcanic eruptions. Other studies, such as Driscoll et al. (2012) and Stenchikov et al. (2006), also find no consistent warming in the Northern Hemisphere or strengthening of the polar vortex associated with winter warming.

To better understand why a strengthening of the polar vortex circulation occurs, several studies have proposed mechanisms that link volcanic eruptions with changes in atmospheric circulation (Robock and Mao, 1995; Robock, 2000; Stenchikov et al., 2002). Despite proposed mechanisms, however, some studies suggest that the prevalence of this response may depend on aerosol forcing (Toohey et al., 2014) or may be insignificant in comparison with the range of natural variability in climate (Polvani et al., 2019).

3 Model description and experimental setup

To investigate the Pinatubo response under different initial conditions we run a large ensemble of simulations with initial conditions sampled using two different sampling schemes. Both sets of simulations are run in accordance with the VolMIP protocol with a pre-industrial atmosphere in GISS model E2.1.

3.1 The model

All model simulations are run in GISS model E2.1 (Kelley et al., 2020) (E2-1-G in CMIP6 archive): a climate model with fully coupled ocean–atmosphere dynamics and in correspondence with CMIP6 protocols. GISS model E2.1 has a horizontal resolution of 2^∘ latitude by 2.5^∘ longitude and 40 vertical layers (which are more densely layered close to the surface and get progressively coarser going upwards into the stratosphere). All ensembles are run with a fully dynamic mass-converting free surface Russel ocean model (Russell et al., 1995) now referred to as the GISS ocean, denoted “G” in GISS E2.1-G. The atmosphere is represented with non-interactive (NINT) aerosols. Thus, ozone and other aerosols are pre-determined by CMIP6- specified model inputs.

The representation of ENSO in GISS E2.1-G for CMIP6 has improved significantly upon E2 (CMIP5) (Schmidt et al., 2014) on correlated global changes in temperature for all ocean representations (including the GISS ocean used here). The model shows a spectral density of ENSO events peaking at a 5-year period (Kelley et al., 2020) showing a slight bias in frequency, although the relative strength of the ENSO cycle is reasonable. E2.1 also shows a higher than average standard deviation in NAO patterns when compared with observations and other models (Orbe et al., 2020) translating to larger variability in the NAO.

3.2 Model simulations and sampling

After a 6000-year control run spin-up, 400 years are sampled to select ensemble years based on initial ENSO and NAO conditions. Prior to sampling, control run simulations are run with a pre-industrial atmosphere for a total of 400 years. Simulations begin in June where the start years are determined by selecting years with specific ENSO and NAO tendencies at the time of peak forcing defined as the average conditions across December (start year), January (start year + 1), and February (start year + 1) as specified by the sampling protocol described in Zanchettin et al. (2016). More details about the sampling protocol and indices for both ENSO and NAO conditions are described by Zanchettin et al. (2016), where we use the Niño 3.4 index to represent ENSO states and a two-box 500 mb geopotential height index to represent NAO (Stephenson et al., 2006). Sampling for the VolMIP simulations is with the same sampling protocol as for other models from the multi-model ensemble shown in Zanchettin et al. (2022). While it is important to note that the sampled states are at the time of peak aerosol forcing rather than at the start of the simulation, we hereafter refer to these climate conditions as initial conditions because they represent the expected climate state in the absence of volcanic forcing. In total, we sample 9 years from each co-condition (a combined cross of three conditions of $\text{ENSO}+/\mathrm{0}/-$ and three conditions of $\text{NAO}+/\mathrm{0}/-$) for a total of 81 VolMIP sampled simulations; each ENSO (and NAO) state thus has 27 members.

In addition to the VolMIP runs sampled from initial conditions, we also sample 50 additional runs randomly from the same control run (henceforth referred to as “random samples”). From these 50 randomly sampled years, 10 overlap with already sampled VolMIP simulation years. Thus, 40 additional simulations (identical to VolMIP simulations except for random sampling) were also run with NINT (non-interactive, specified aerosols, and atmospheric chemistry) atmosphere and the GISS ocean for a total of 121 simulations. Initial conditions in the control run are approximately randomly distributed and show some skew in ENSO states, but no correlation between ENSO and NAO states (Fig. 1).

Figure 1 Initial climate conditions for 400 years in the GISS control model run. Each year in the 400-year control is plotted with ENSO (Niño 3.4) index on the x axis and NAO (z-500) index on the y axis. VolMIP sampled simulations are denoted by red circles, additional randomly sampled years are blue (+), and all other years in the control run are black dots. Histograms show a distribution of initial conditions within the model control run: ENSO conditions exhibit a positive skew consistent with model biases addressed by Kelley et al. (2020).

After ensemble years are selected according to the two sampling schemes (VolMIP and random samples), volcanic simulations are run in GISS-E2-1-G in accordance with VolMIP protocol (Zanchettin et al., 2016). Volcanic aerosols are prescribed based on CMIP6 Pinatubo aerosol climatology (Thomason et al., 2018) as a function of height, latitude, and time beginning on the sixth month of the simulation (June) to emulate the 1991 Pinatubo eruption.

3.3 Data analysis and anomalies

We examine how the climate response under volcanic conditions differs from control conditions under two different widely used methods to represent the climate anomaly of this type of short-lived climate perturbation. First, we process results by computing anomalies from the equivalent control period (response = perturbed − control), hereafter referred to as “paired anomalies” as defined in Zanchettin et al. (2022). With paired anomalies each ensemble is analyzed as an anomaly from control conditions with the same initial climate condition, thereby excluding the effect of ongoing unperturbed climate variability, initial ENSO, and NAO state (Pausata et al., 2015). The second approach we examine calculates anomalies from a climatological control period (hereafter “climate anomalies”) which instead measure deviation from the control run and thus include ongoing climate variability. Specifically, we use two definitions of climatological anomalies to best contextualize our results with other studies. The first climatological anomaly takes the difference between each ensemble member and the full control including seasonality as defined by Zanchettin et al. (2022). Second, for comparison of Northern Hemisphere anomalies we additionally calculate a second climate anomaly by subtracting the average condition of the 5 years prior to the selected eruption year as used in other studies examining the climate response to volcanic forcing (Polvani et al., 2019; Stenchikov et al., 2006; Driscoll et al., 2012). However, in order to discuss responses in the absence of background conditions, we present predominately paired anomalies except where climatological anomalies produce significantly different results.

For analysis of VolMIP sampled simulations (see Sect. 1.2), we compute means not only with all 81 VolMIP ensembles but also with subsets of ensembles grouped according to their initial ENSO and NAO ensembles. For example, we look at the ensemble mean of all ensembles that were sampled with positive ENSO initial conditions (+ENSO) and compare these 27 ensemble members to the 27 sampled ensembles which show neural (0 ENSO) and negative (−ENSO) conditions. We follow this same nomenclature for analysis of differences grouped by initial NAO conditions, and thus when discussing the response of a given group (e.g., −ENSO) we refer to the mean paired anomalous response of that ensemble group unless otherwise noted.

Where applicable, we also test the statistical significance of differences between VolMIP ensemble groups of different ENSO and NAO initial conditions. In these cases, we compare the values from the 27 ensemble members using an analysis of variance (ANOVA) test and report the corresponding p value to represent the statistical difference ensemble groups. For the case of displaying spatial differences in the surface temperature response, we perform a Student's t test between positive and negative NAO ensemble groups, presenting the mean surface temperature response only at model grid cells where the p value from the t test between positive NAO and negative NAO groups is <0.05.

4 Results

For all ensemble members, we find that radiative forcing impacts and surface temperature impacts are consistent with previous observations and modeled studies of 1991 Pinatubo-sized volcanic eruptions (Schmidt et al., 2018), with ensemble mean forcing peaking at −3.27 W m⁻² the December after the eruption. In comparison with other models in VolMIP, we note that GISS E 2.1 does display a faster increase of radiative anomalies (Zanchettin et al., 2022). However, between our different ensemble members, there is a little variation in the evolution of the radiative response to the prescribed volcanic forcing (see Fig. S1 in the Supplement).

Changes in radiative forcing are also accompanied by a reduction in global surface temperature peaking at −0.35 ^∘C the first spring after the eruption. Further analysis of surface temperature anomalies both globally and in the tropics are presented in Fig. S2.

4.1 ENSO Response

Figure 2 shows the monthly Niño 3.4 index (filtered to remove the seasonal signal) for positive, negative, and neutral VolMIP ensembles with paired anomalies. Positive, neutral, and negative ENSO ensemble groups all show negative, La Niña-like sea surface temperature anomalies in the first post-eruptive winter. Negative sea surface temperature anomalies, however, are strongest for +ENSO and 0 ENSO conditions with mean peak decreases of −0.61 and −0.67 ^∘C, respectively, consistent with changes in tropical surface temperature. In +ENSO ensembles the sea surface temperature relaxes toward mean temperatures from warmer than average conditions. −ENSO ensembles show little variation between control and perturbed simulations, suggesting that a cooler-than-average tropical sea surface temperature will be affected little by volcanic perturbations. The lack of a clear ENSO response for La Niña ensembles is consistent with Predybaylo et al. (2017) where La Niña ensembles showed no significant ENSO anomaly. Unlike results of earlier studies (Pausata et al., 2020; Khodri et al., 2017; Predybaylo et al., 2017), we find no El Niño anomalies in these simulations. Our findings do, however, support the idea that ENSO response is dependent on pre-conditioning or initial conditions in the tropical Pacific (McGregor et al., 2020).

We do not differentiate here between central Pacific and eastern Pacific El Niño events as in Predybaylo et al. (2017). The small inter-ensemble spread for +ENSO simulations, however, suggests that there is little difference between the two in our model representation. We additionally perform an equivalent analysis of the relative sea surface temperature (RSST) Niño 3.4 index as proposed by Khodri et al. (2017) and presented in other studies. The RSST Niño 3.4 anomalies (Fig. S4) also show consistent cooling in the first post-eruptive winter with little difference in the strength of cooling between different initial ENSO conditions. This suggests that in the GISS model the Niño 3.4 region experiences elevated levels of cooling relative to cooling throughout the tropical region. The RSST Niño 3.4 anomalies on average also show a slight warm anomaly in the second post-eruptive year, as noted in Zanchettin et al. (2022) where a subset of these ensembles was analyzed with other VolMIP models.

Figure 2 Monthly, seasonally de-trended time series of changes in the monthly Niño 3.4 ENSO index under different ENSO initial conditions and for the full VolMIP ensemble. Positive and negative ensembles show a relaxation of the index response toward mean climatological sea surface temperature conditions. Red shading shows the 95% confidence interval for the anomalous response from control conditions.

4.2 Northern Hemisphere response

We now turn to the response in the Northern Hemisphere, which has also been widely discussed between model studies. Specifically, we consider how initial NAO conditions impact the Northern Hemisphere climate response, particularly in the first winter.

4.2.1 NAO response

Figure 3 shows the monthly NAO index (based on 500 mb geopotential height) throughout the 3-year simulation period for positive, neutral, and negative NAO groups. Regardless of the initial condition of the NAO, the NAO relaxes toward mean conditions. For +NAO ensembles, this is shown by a decrease in geopotential height, peaking at 69.7 mb in the February after the eruption. Likewise, −NAO ensembles show increased geopotential height toward mean conditions by an average of 88.5 mb peaking the first February after the eruption. Neutral NAO (0 NAO) ensembles have no statistically significant response, with the confidence interval of these ensembles showing significant variability in mean geopotential height. When looking at groups of specific initial NAO conditions, it is clear that for +NAO and −NAO conditions the geopotential height is relaxed toward mean conditions. When looking at the full 81-member ensemble (Fig. 3), however, there is a tendency toward positive NAO anomalies in the first post-eruption winter in the model as noted by Zanchettin et al. (2022).

These findings suggest that in our simulations, regardless of initial conditions, there is a relaxation of extreme positive or negative NAO conditions otherwise present in control runs. For negative NAO ensemble years, this causes an anomalous strengthening of the pressure dipole between the Azores high and Icelandic low regions when eruptions occur under negative NAO conditions. The opposite is true for +NAO simulations, where the dipole between these two pressure systems appears to weaken in comparison with control conditions the first winter after.

Given this robust change in pressure in the North Atlantic during the first winter, we continue to discuss how these changes are seen through other polar dynamic pathways.

Figure 3 Monthly, seasonally de-trended time series of changes in the monthly NAO index under different NAO initial conditions. Red shading shows the 95 % confidence interval for the anomalous response from control conditions. +NAO ensembles show a robust decrease in the NAO index in the first winter ($t=\mathrm{12}-\mathrm{14}$), while negative NAO ensembles show a robust increase.

Figure 4 Regressions of the (a) Equator-to-pole temperature gradient (temperature at 50 mb at the Equator; temperature at 50 mb at the pole) vs. stratospheric polar vortex strength (u₅₀ at 61^∘ N) and (b) winter warming vs. polar vortex strength, all in the first post-eruptive winter. All 81 VolMIP ensemble members are plotted with shape and color corresponding to the initial NAO phase. R² values are displayed for each regression.

4.2.2 Polar dynamic changes

Modeled changes in the North Atlantic geopotential height dipole (quantified by the NAO index) are accompanied by other changes in zonal winds and atmospheric temperature. Specifically, we find that the polar vortex strength (defined as the zonal mean wind at 10 hPa and 60^∘ N as in Polvani et al., 2019) show positive anomalies from control conditions for −NAO ensembles and negative anomalies for +NAO ensembles (Fig. 4). These variations in polar vortex strength between NAO groups are statistically significant with a p value of $\mathrm{2.15}\times {\mathrm{10}}^{-\mathrm{12}}$.

These changes in zonal winds are also accompanied changing patterns in 500 mb geopotential height (∼5.5 km) over the polar region (60–90^∘ N). For negative NAO ensembles, geopotential height decreases by an ensemble average of 100 mb (∼15.5 km) in the first winter. There are also moderate increases in geopotential height in the Atlantic basin near the Azores (20–55^∘ N, 90^∘ W–60^∘ E), averaging around 50 mb. These changes in pressure and wind are indicative of strengthening of the polar vortex and a positive phase of the NAM. The opposite occurs for +NAO ensembles, consistent with a decrease in the strength of the polar vortex and a negative phase of the NAM. Mean changes in geopotential height for each NAO group are presented in Fig. S5.

4.2.3 Atmospheric temperature

Analysis of changes in temperature in the lower stratosphere illustrate how volcanic eruptions impact atmospheric temperature under different initial conditions. Specifically, we examine temperature anomalies using the modeled microwave sounding unit temperature metrics in the lower stratosphere (MSU TLS). The MSU temperature metric is commonly used as a remotely sensed temperature data metric based on height; however, here we present an equivalent modeled metric in E2.1 (Miller et al., 2020). Figure 5 shows the anomaly in MSU temperature in the lower stratosphere for the first boreal winter across latitude and time. All ensembles exhibit robust tropical tropospheric warming peaking at 2.5 ^∘C, tapering off toward the south pole.

The temperature anomaly north of 60^∘ N, however, varies significantly between simulations. +NAO ensembles show an anomalous warming in the stratosphere reaching an average of 2.5^∘ at the North Pole. −NAO ensembles show stratospheric cooling anomalies reaching an average of 5.8 ^∘C at the North Pole. Neutral NAO ensembles again fall between these extremes falling close to the mean.

For all ensembles, volcanic forcing smooths out meridional temperature gradients in the first winter that are present in control conditions. Thus, −NAO ensemble simulations, which would normally have a weak meridional temperature gradient, increase the high northern latitude gradient in the first winter after the volcanic eruption. The opposite occurs for +NAO ensembles, where higher than average temperature gradients are decreased to mean conditions (Fig. S3).

Figure 5 Panel (a) shows the change in lower stratospheric temperature and (b) shows the unperturbed lower stratospheric temperature derived from microwave sounding units in the first winter (December–February) after the eruption. Shading denotes the 95 % confidence interval for each ensemble group. Tropical stratospheric warming occurs for all ensembles; however, the high northern latitude response varies greatly between different initial conditions. −NAO ensembles show cooling in the high latitude lower stratosphere while +NAO ensembles show significant warming.

Changes in the polar stratospheric temperature drive the strength of the Equator-to-pole temperature difference. −NAO ensembles drive an increase in the Equator-to-pole temperature difference driven both by warming of the equatorial lower stratosphere and cooling of the high latitude lower stratosphere. +NAO ensembles, on the other hand show little change in the Equator-to-pole temperature difference as both the equatorial and polar lower stratosphere experience warming. Neutral NAO ensembles fall somewhere in the middle with a moderate increase in the Equator-to-pole temperature difference. To further investigate if an enhanced Equator-to-pole temperature difference correlates with an increased polar vortex circulation, we use a simple regression as done by Polvani et al. (2019) with each of our 81 VolMIP ensembles. Figure 4 shows that changes in the Equator-to-pole temperature gradient in the first winter strongly correlates with increased polar vortex strength in the first winter. There is also a correlation between the observed winter warming anomaly and polar vortex strength (R²=0.40) indicating that a strengthening of the polar vortex often corresponds with winter warming. The correlation between vortex strength and winter warming is stronger than in Polvani et al. (2019), which could suggest that larger ensembles are required to find a significant signal but does not suggest that a strengthened polar vortex alone is the cause of observed winter warming. We also note that there is also a significant difference in the Equator-to-pole temperature gradient between NAO groups with a p value of $\mathrm{9.68}\times {\mathrm{10}}^{-\mathrm{8}}$.

4.2.4 Winter warming

Having discussed dynamic changes in the Northern Hemisphere, we now discuss the strength of the winter warming response across different initial conditions. Figure 6 shows the mean surface temperature anomaly in the first boreal winter (December–February) after the eruption for positive, negative, and neutral ensemble groups. Most areas (where shading is gray) do not experience any robust difference between initial NAO phases. Northern Eurasia and Greenland, however, have significantly different responses between $+/-\text{NAO}$ conditions. +NAO ensembles experience cooling over Eurasia and warming over Greenland. −NAO ensembles show the opposite, with significant warming over Eurasia and cooling over Greenland. Neutral NAO ensembles show a weak warming signature similar to that of −NAO ensembles. The winter warming signature (measured as a mean of 27 ensembles) is strong only for the −NAO group; however, we note that both +NAO and −NAO ensembles under perturbed conditions have a trend toward the surface temperature displayed during 0 NAO conditions.

Figure 6 Average surface temperature paired anomaly in the first (December–February) winter after the eruption for positive (a), neutral (b), and negative (c) NAO ensembles. Gray areas for positive and negative ensembles denote confidence below 95 % in the difference between positive and negative ensemble groups.

While −NAO ensemble group means show a significant winter warming response, now we look at variation within NAO groups. Figure 7 shows a box plot for the winter warming in Eurasia (40–70^∘ N, 0–150^∘ W) of simulations grouped by NAO phase, with all VolMIP ensembles, and for the 50 randomly sampled simulations for comparison. We also include anomalies from historical conditions for direct comparison with other studies (Polvani et al., 2019; Driscoll et al., 2012).

When using paired anomalies −NAO ensembles show mean warming over Eurasia, with few ensembles showing a negative anomaly. +NAO ensembles all experience a cooling temperature anomaly in Eurasia. The neutral group of ensembles has a mean around 0^∘ of warming but is slightly skewed to a positive temperature anomaly. An ANOVA test shows there is a statistically significant value between paired anomaly ensemble groups with different initial NAO conditions in the VolMIP ensembles with an F statistic value of 22.78 and a p value of $\mathrm{1.62}\times {\mathrm{10}}^{-\mathrm{8}}$. Plotting all VolMIP ensembles together shows a large variation in the temperature response due to including all initial conditions together. The randomly sampled runs have a distribution similar to the neutral NAO ensemble groups, suggesting that extreme initial conditions, such as very negative NAO or positive NAO phases, are less common in the climate system than in our sample.

Climatological anomalies show no significant forced response, contrary to the paired anomalies for +NAO and −NAO groups. This suggests that paired anomalies are influenced by the sampled conditions in the unperturbed control. These sampled states of NAO are evident as paired anomalies show cooler than average conditions for −NAO ensembles and warmer than average conditions for +NAO ensembles. There is, however, no significant difference between the perturbed (with volcanic forcing) and control (with no volcanic forcing) winter warming response for all ensemble members (All VolMIP and Random Samples) or for neutral NAO ensemble members.

In addition to decreasing variability in the response, ensemble groupings also impact the probability of observing warming in the model (when considering anomalies taken from control conditions). Table 1 shows the probability of simulations showing winter warming (calculated using paired anomalies) for varying initial conditions calculated using the 81 VolMIP runs and the 27-member initial condition groupings. While not all −NAO ensembles show a winter warming response, the probability of observing winter warming increases greatly for negative NAO initial conditions in comparison with neutral and negative conditions when using paired anomalies. The probability for observing a significant winter warming response in any of the 81 VolMIP samples is low (32 %). The probability of warming given −NAO initial conditions, however, is higher (60 %). Thus, while these simulations still show the large variation in surface temperature responses the first winter after the eruption, initial conditions impact how likely a paired anomalous winter warming response is in a large group of ensembles.

Table 1 Winter warming probabilities (computed with VolMIP-sampled simulations).

5 Summary and discussion

Initial ENSO conditions lead to significant differences in the Niño 3.4 sea surface temperature anomaly through the first boreal winter into the spring after the eruption. The temperature decrease is present in the ensemble mean but strongest for ensembles with positive and neutral ENSO initial conditions. In general, we find no signature of an enhanced El Niño-like anomaly in the first winter, as has been suggested in other studies (Pausata et al., 2020; Khodri et al., 2017; Predybaylo et al., 2017); however, there is a slight warm anomaly following the initial cooling after the eruption.

The response of the Northern Hemisphere varies significantly between ensembles with different initial NAO conditions both in sign and strength of responses. Winter warming anomalies occur with increased probability for ensembles with −NAO initial conditions with 60 % of ensembles showing a warming response in the first winter. This warming response corresponds with an anomalous decrease in polar lower stratospheric temperature in the first winter for −NAO ensembles, causing an increased temperature gradient between the Equator and poles. A simple regression analysis shows that positive temperature gradient anomalies are correlated with an increased strength of the stratospheric polar vortex. −NAO ensembles also exhibit decreased geopotential heights and increased westerly zonal wind circulation that are consistent with this strengthening polar vortex following the eruption. There is a weak correlation between strengthening of the stratospheric polar vortex and the winter warming response in the first winter, although this correlation does not suggest the lack of other response pathways. In general, +NAO ensembles show the opposite anomalous patterns from control conditions and neutral ensembles show some weak warming and vortex strengthening anomalies. Thus, while a −NAO phase does not guarantee winter warming resulting from the Equator-to-pole temperature difference, it does highly increase the probability that winter warming will occur. These polar dynamic changes also coincide with a smoothing of meridional temperature gradients from control conditions. This response could be model dependent or a result of the specific way that the Pinatubo forcing is prescribed in the simulation with non-interactive aerosols.

Figure 7 Simulated boreal winter (December–February) warming response in Eurasia the first year after the eruption for both control (no volcanic forcing) and perturbed (with volcanic forcing) runs. Control anomalies (gray) are taken from the mean winter surface temperature for the 5 years prior to the eruption. Perturbed anomalies are shown both with anomalies from control conditions (red) and from historical conditions (blue). Box plots are shown for all 81 simulations (All VolMIP), for each initial condition group (Negative, Neutral, and Positive), each with 27 simulations and for all 50 randomly sampled runs.

For all simulations, the monthly NAO index in the first winter relaxes toward mean conditions. This means that for both +NAO and −NAO ensembles, there is a sudden anomalous change in pressure in the North Atlantic after the eruption. Thus, the anomalous strengthening of the polar vortex from control conditions could be due to the sudden relaxation of the NAO anomaly in the first winter relative to initial conditions. The strengthening of the stratospheric polar vortex resulting in winter warming thus only occurs when the model would have otherwise experienced weak vortex circulation in the absence of volcanic forcing. The anomalous response is significantly impacted by these extreme initial conditions that were sampled from our control conditions. To test this impact, we also compare our 27-member $\text{NAO}+/\mathrm{0}/-$ ensemble groups to a 50-member randomly sampled ensemble group. The randomly sampled ensemble group shows anomalies most similar to the neutral NAO ensemble group, suggesting that strong anomalies due to extreme NAO initial conditions are less common in a representative sample. While extremely negative phases are most likely to experience winter warming such extremes are less common in the real world, possibly explaining why warming is only sometimes observed in model simulation ensembles. These extremes in initial conditions can contribute significantly to ensemble variation, particularly with a small amount of ensemble members or when ensembles are sampled with a bias in initial conditions.

We tested the impact of using paired anomalies by comparing our analysis with climatological anomalies used in other studies (Polvani et al., 2019; Driscoll et al., 2012). These climatological anomalies, which take reference from mean climate conditions, show no statistically significant forced response for our ensemble members. The difference in responses between anomalies demonstrates how the choice of anomalies can significantly impact the modeled response, particularly when extreme initial conditions are present as initial climate conditions can significantly influence the strength of a given response. For example, when analyzing the winter warming response under varying NAO conditions, ensembles with a strong −NAO condition in the control run relaxed toward mean NAO conditions under perturbed volcanic aerosol runs. Thus, when using paired anomalies the strength of the winter warming response is biased due to lower mean air temperatures in the control simulation. When using historical anomalies, we see no significant warming response for the same perturbed runs, as air temperatures are typical of the historical-mean climate state (neutral NAO conditions). The difference between the modeled responses under paired and historical anomalies was also highlighted by Zanchettin et al. (2022), where the choice of anomaly was shown to impact some ensemble mean responses.

Here simulations have been constrained to examine the climate response with a protocol that eliminates some sources of variability. In particular, VolMIP compliant simulations used here are run with non-interactive (NINT) aerosols and represent pre-industrial conditions. Thus, they cannot be directly compared to Pinatubo simulations which have industrial greenhouse gases and other constituents in the atmosphere. These runs also do not account for changes in ozone concentration observed after eruptions which may also influence changes in stratospheric circulation (Stenchikov et al., 2002). The NINT atmospheric representation also dictates that aerosols evolve exactly as prescribed, making aerosols insensitive to states of the stratosphere and troposphere. Other runs with interactive aerosols are necessary to understand if the dynamics of these responses are dependent on the specific prescription of volcanic aerosols.

Furthermore, the current GISS model E2.1 does not have a realistic representation of some key atmospheric components such as the quasi-biennial oscillation (QBO) (Rind et al., 2020) that could also play a role in the observed circulation responses (Stenchikov et al., 2004). Changes in the QBO could also influence the strength of the polar vortex circulation as easterly phases (such as those during the Pinatubo eruption) are likely to cause a decrease in the stratospheric polar circulation (Holton and Tan, 1980). Here we have used the GISS E2.1-G CMIP6 compliant runs; however, in a future study, we hope to examine this response in the new GISS model E2.2+ which has a higher vertical resolution and a better representation of the QBO.

Overall, we find that initial ENSO conditions have a small effect on surface temperature and ENSO response as a cooling, La Niña-like, anomaly in the tropical Pacific occurs for each ensemble in the first post-eruptive winter. The initial state of the NAO, on the other hand, varies the anomalous response by relaxing initial conditions in the first winter to a neutral NAO phase. If a volcanic eruption occurs during a normally −NAO phase, these changes in turn increase the probability of observing a winter warming response from control conditions in the first post-eruptive winter. For extremes in initial NAO conditions, changes in the Northern Hemisphere are the most robust. While often these extremes are uncommon, they likely contribute to inter-ensemble variation and thus uncertainty in predicting the climate's response to volcanic eruptions. When the forced winter warming is defined as the average of a large ensemble (including all initial conditions), however, the response is insignificant (mean with an ensemble spread around zero). The prevalence and strength of this anomaly is influenced both by extremes in initial conditions and how anomalies are taken (either paired or climatological anomalies).

This study highlights the takeaway from Zanchettin et al. (2022) that initial conditions have a significant impact on post-eruptive anomalies. Specifically, our expanded 81-member ensemble shows that in GISS E2.1-G, there is a statistically significant difference in how the Northern Hemisphere responds to volcanic forcing under different initial conditions of NAO when using paired anomalies. This implies that for future studies using small ensembles, responses such as winter warming may be significantly impacted by extremes in initial NAO conditions, and climatological anomalies should also be explored to show robustness in modeled responses.

6 Conclusions

The climate response to large, Pinatubo-type volcanic eruptions is variable between models and has been discussed here in relation to GISS model E2.1-G. We focus on two responses which have been studied both with observational and modeling studies: the ENSO response and Northern Hemisphere response in the first winter. A total of 121 ensembles were run in the GISS E2.1-G model to examine how initial ENSO and NAO conditions impact the modeled climate response. Our experimental setup uses a pre-industrial model with prescribed aerosols and took anomalies from an equivalent control period run rather than an historical climate period, allowing us to filter out initial climate variability and look only at the response due to volcanic sulfate aerosols. We find that ensembles with different initial NAO conditions have significantly different anomalous climate responses in the first Northern Hemisphere winter. In particular, years which would be in +NAO or −NAO conditions are relaxed to mean NAO conditions under volcanic forcing. This creates an anomalous negative and positive winter warming response for +NAO and −NAO ensembles, respectively. Ensembles with different initial ENSO conditions, however, show similar anomalies between different initial phases. Thus, inter-ensemble variation caused by initial conditions is significant particularly when looking at the first Northern Hemisphere winter response.