1 Introduction

Volcano eruptions are a dominant external forcing for triggering interannual climate responses through injecting large amounts of sulfur dioxide into the stratosphere and perturbing the radiative balance (Robock 2000; Hegerl et al. 2003; Schneider et al. 2009; Timmreck 2012; Myhre et al. 2014; Schurer et al. 2014). In the instrumental observations for the period of 1870–2010, an El Niño occurred in the first winter following the 1902 Santa Maria, 1963 Agung, 1982 El Chichón, and 1991 Pinatubo eruptions, but not after the 1883 Krakatau eruption (Khodri et al. 2017). These observations indicate that explosive eruptions may have important impacts on initialization and development of an El Niño, which is the positive phase of El Niño-Southern Oscillation (ENSO). This hypothesis, however, cannot be convincingly confirmed by instrumental observations due to the small sample size (Graf et al. 2014).

Fortunately, some long-term reconstructions have been developed; and we can use them to test the hypothesis with more samples. Previous studies suggested that large tropical eruptions increase the likelihood of an El Niño within two years after the eruption (Handler 1984; Adams et al. 2003; Emile-Geay et al. 2008; D’Arrigo et al. 2009; McGregor et al. 2010; Wilson et al. 2010; Li et al. 2013). Recent analyses based on an 1100-year reconstruction suggested that an El Niño mainly occurs in the first winter after a large tropical eruption (Liu et al. 2018a, b).

Many state-of-the-art models can successfully simulate this El-Niño-like response to large tropical eruptions (Mann et al. 2005; McGregor et al. 2010; Wilson et al. 2010; Ohba et al. 2013; Maher et al. 2015; Lim et al. 2016; Khodri et al. 2017; Predybaylo et al. 2017; Stevenson et al. 2017; Liu et al. 2018a; Pausata et al. 2020) and to Northern-Hemispheric eruptions (Pausata et al. 2015, 2016; Stevenson et al. 2016; Liu et al. 2018a), but some models cannot (McGregor and Timmermann 2011; Ding et al. 2014; Maher et al. 2015; Wang et al. 2017; Liu et al. 2020). The physical linkage between volcano eruptions and El Niño development has been studied in the past. After a large tropical eruption, an El Niño can be excited by the ocean dynamical thermostat mechanism (Clement et al. 1996) under the overall tropical cooling in some models (Mann et al. 2005; Maher et al. 2015). An El Niño can also be induced by the equatorward migration of the Intertropical Convergence Zone (ITCZ) due to asymmetrically distributed volcanic-induced cooling (Pausata et al. 2015, 2016, 2020; Stevenson et al. 2016) or to less evaporation in the subtropical cloud-free region caused by a cooling associated with the reduction of net surface shortwave radiation (Lim et al. 2016).

Excitation of a westerly anomaly over the western-to-central equatorial Pacific (WCEP) after a tropical eruption is also found to be an important trigger for El Niño development (Ohba et al. 2013; Khodri et al. 2017; Predybaylo et al. 2017; Liu et al. 2018a), mainly through the Bjerknes feedback (Bjerknes 1969). An initial westerly anomaly is suggested to be excited by a land-sea thermal contrast around the Maritime Continent (MC) (Ohba et al. 2013), suppressing tropical precipitation (Liu et al. 2018a), or by a Kelvin wave response to a weakened West African monsoon (Khodri et al. 2017). However, these findings were obtained based on analyses of individual models. It is still unclear whether a WCEP westerly anomaly is a common response to tropical eruptions in models or not.

This study examines tropical responses in multiple models, focusing on westerly anomalies excited by tropical eruptions. It is not clear whether the WCEP westerly anomaly guarantees the occurrence of an El Niño in these models; thus, we further study the simulation of the El Niño after a WCEP westerly anomaly excitation. These results should help enhance our understanding of climate system responses to one of the most important natural forcings, and improve climate models in the future.

2 Data, model, and methods

2.1 PMIP3 simulations

In this study, we analyze the response to large tropical eruptions from all transient climate simulations of the last millennium in the Paleoclimate Modeling Inter-comparison Project Phase 3 (PMIP3) archived from 850 to 1850 AD (Braconnot et al. 2011, 2012). The 10 PMIP3 models used here are listed in Table 1. More details on these models can be found in Masson-Delmotte et al. (2013). Different forcing datasets are used in these full-forcing millennium simulations, following the protocols of the PMIP3 (https://pmip3.lsce.ipsl.fr/wiki/doku.php/pmip3:design:lm:final), as discussed in Schmidt et al. (2011, 2012). To force these millennium simulations, two sets of reconstructions for stratospheric volcanic aerosols were derived from either Crowley et al. (2008), which is based on 13 Greenland and Antarctic ice-core records, and from Gao et al. (2008), which is derived from 32 Arctic and 22 Antarctic ice-core records (Table 1). Only the IPSL-CM5A-LR uses the volcanic forcing of Ammann et al. (2007), which was established by converting ice core aerosol proxies to latitudinal and temporally varying atmospheric aerosol fields (Atwood et al. 2016).

Table 1 Description of PMIP3 and CESM models, volcanic forcing datasets, number of runs, period, and references
Full size table

2.2 CESM simulations

In this study, we also use the outputs from 850 to 1850 AD of the “full forcing” and “volcanic only” ensembles in the Community Earth System Model-Last Millennium Ensemble (CESM-LME) (Otto-Bliesner et al. 2016), which consist of 10 and 5 ensembles, respectively (

Table 1). In “volcanic only” ensembles, only volcanic forcing perturbations are imposed, while in the “full forcing” ensembles all external forcing perturbations, such as solar variability, volcanic eruptions, and greenhouse gases (GHGs) forcings, are imposed. The volcanic forcing field for both sets of ensembles is the monthly, zonally averaged, meridional distributions of stratospheric aerosols from Gao et al. (2008).

We also use the CESM-Large Ensemble (LE) experiment with 41 ensemble members from 1920 to 2005 (Kay et al. 2015) to investigate model responses to the 1963 Agung, 1982 El Chichón, and 1991 Pinatubo eruptions. All the CESM-LEs have the same external forcings following the fifth phase of the Coupled Model Intercomparison Project (CMIP5) design protocol (Taylor et al. 2012), which is forced by time-evolving observed GHGs, aerosols, and other radiative forcing from 1920 to 2005 (Lamarque et al. 2010). Volcanic forcing is incorporated in this experiment using the dataset of Ammann et al. (2003). In the CESM-LE, the strength and the volcanic aerosol distribution are idealized and symmetric around the equator for all eruptions considered in this study. Details of this CESM-LE experiment can be found in Kay et al. (2015).

2.3 Observational and proxy data

To analyze the observations of the three recent large eruptions of 1963 Agung, 1982 El Chichón, and 1991 Pinatubo (Sato et al. 1993), instrumental sea surface temperature (SST) dataset are obtained from the Hadley Centre Ice and SST version 1 (HadISST1) (Rayner et al. 2003) and un-interpolated version 2 (HadISST2) (Rayner et al. 2006). The dataset of the Extended Reconstructed SST dataset version 5 (ERSSTv5) (Huang et al. 2017) is used as a reference. Monthly-mean 850 hPa wind is obtained from the National Centers for Environmental Prediction (NCEP)/National Center for the Atmospheric Research (NCAR) Reanalysis 1 (NCEP-R1) (Kalnay et al. 1996). The Niño 3.4 index is defined as the area averaged SST over (170°–120° W, 5° S–5° N).

We select the 10 largest tropical eruptions with different strengths and durations from 850 to 1850 AD. These tropical eruptions as well as three recent large tropical eruptions during the last century all had larger aerosol column density or aerosol optical depth (AOD) in the tropics than in the mid-to-high latitudes (Figs. S1 and S2) and lasted for 2–3 years (Fig. S2). The largest injection occurred in 1258 (Oppenheimer 2003), and was associated with the eruption of the tropical Samalas volcano, part of the Mount Rinjani Volcanic Complex, Indonesia (Lavigne et al. 2013). This event was followed by three strong tropical eruptions in 1268, 1275 and 1284. These eruptions were suggested to have triggered the Little Ice Age (Miller et al. 2012; Slawinska and Robock 2018). For each individual simulation, the actual volcanic forcing dataset of Crowley et al. (2008) or Gao et al. (2008) is used to identify the eruption dates. In Gao et al. (2008), the volcanic forcing is specified to peak in June. In Crowley et al. (2008), three of these 10 largest eruptions, i.e., the 971, 1809, and 1835 eruptions (Table 2), are specified to peak in January; and the others, in June. Here, the eruption year 0 is defined as the peak loading year. Details of the sulfate injections of these 10 eruptions are shown in Table 2, including the actual dates that appear in the two different forcing time series used for each individual simulation. The eruption season, which has been found to affect atmospheric circulation responses in simulations (McGregor and Timmermann 2011; Predybaylo et al. 2017; Stevenson et al. 2017; Sun et al. 2019), is not considered in this work, which should be studied in the future.

Table 2 Ten largest tropical volcanic eruptions during 850–1850 AD
Full size table

2.4 Model description and experiment design

Under the volcanic forcing, land becomes cool much faster than ocean. To study the role of land cooling in exciting the tropical circulation, simulations are performed using the CESM 1.2.2 (Hurrell et al. 2013). The CESM configuration we used is composed of Community Atmosphere Model version 5.0 (Neale et al. 2010) and Community Land Model version 4.5 (Oleson et al. 2013) components, with realistic monthly SST and sea ice specified from 1979 to 2013, which means that we follow the Atmospheric Model Intercomparison Project (AMIP) protocols, and only the atmosphere and land model are integrated dynamically. The model resolution is 0.9° (latitude) × 1.25° (longitude) and 26 vertical levels.

Based on the selected 10-year control run from 1980 to 1989, we conduct seven sets of sensitivity experiments by doubling the land albedo over the globe, extratropical (0–360°, 10°–90° S, 30°–90° N), tropical (0–360°, 10° S–30° N), South–South East Asia (SSEA, 70–154° E, 10° S–30° N), Asian monsoon (AM, 70–130° E, 10°–30° N), MC (96–154° E, 10° S–10° N), and tropical African (0–52° E, 10° S–30° N) regions to reduce irradiance (Fig. S3). Here, 10° S–30° N is used for dynamic tropical region from June to November. The change of albedo can be achieved by perturbing the surface upward scattered flux, and more details can be found in Oleson et al. (2013). Each set composites 10 integrations and each one-year long. We generate the individual ensemble members by perturbing the initial conditions of the specific year; perturbations are constructed by replacing the state of the atmosphere on January 1st from 1980 to 1989 in the control run. Average of these 10 ensembles is calculated to study the land surface cooling effect.

2.5 Gill model

To study the circulation response to a tropical precipitation anomaly of the first baroclinic mode, the Gill (1980) model is a good tool. The role of the background mean flow is not considered here. Choosing velocity scale of \(C = 50 {\text{ms}}^{ - 1}\) (the speed of the lowest internal gravity wave), length scale of \(\sqrt {C/\beta } = 1500\; {\text{km}}\) (the equatorial Rossby deformation radius), and time scale of \(\sqrt {1/C\beta } = 8.5\; {\text{h}}\), where \(\beta = 2.3 \times 10^{ - 11} {\text{ m}}^{ - 1} {\text{s}}^{ - 1}\) represents the leading-order curvature effect of the Earth at the equator, the non-dimensional Gill model can be written as

$$- yv + \phi_{x} = - \varepsilon u$$
$$yu + \phi_{y} = - \varepsilon v$$
$$(u_{x} + v_{y} ) = - Q - \varepsilon \phi ,$$
(1)

where \(u\) and \(v\) denote lower-tropospheric zonal \(x\) and meridional \(y\) winds, respectively; y is the meridional distance from the equator; and \(\phi\) is geopotential height. \(Q\) is condensational heating associated with a precipitation anomaly, and \(\varepsilon\) is the damping for both momentum and heat. In this work, a damping scale of 3.5 days is used, as in Gill (1980). Sensitivity experiments on this damping scale will also be performed.

2.6 Methods

To understand the volcanic forcing-induced hydrological change, the perturbation energy budget is applied here (Mitchell et al. 1987; Allen and Ingram 2002). For the low-frequency variations with radiative-convective equilibrium, global mean precipitation should be balanced by evaporation. Following the analysis of Wu et al. (2013a), the perturbation energy budget can be written as

$$L\Delta P = \Delta R_{LT} + \Delta R_{ST} + \Delta R_{LS} + \Delta R_{SS} + \Delta Q_{S} ,$$
(2)

where \(L = 2.5 \times 10^{6} {\text{ Jkg}}^{ - 1}\) is the latent heat of condensation and \(P\) is the precipitation change. \(R_{LT}\) and \(R_{ST}\) are the net longwave and shortwave radiation change at the top of the atmosphere (TOA), respectively, and \(R_{LS }\) and \(R_{SS }\) are the net surface longwave and shortwave radiation change, respectively. \(Q_{S}\) is the sensible heat flux change. Δ denotes the global average. In this work, the positive value denotes upward at the TOA and downward at the surface.

The superposed epoch analysis (SEA) (Haurwitz 1981) is used to evaluate the influence of explosive tropical volcanoes on global cooling and related atmosphere–ocean interaction. The influence of internal variability is expected to be removed by performing a SEA composite on the 10 largest tropical eruptions in each model, especially in CESM since it has more samples than the other models. The response to each eruption can also be isolated by performing a SEA composite on the 41 ensemble members of CESM-LE for each eruption. To isolate climate responses to an eruption from the background noise, we obtain the anomaly by removing the average of the climatology and annual cycle based on the five years preceding the eruption. Significance of the SEA is calculated by the bootstrapped resampling method with 10,000 random draws from the study period (Adams et al. 2003). We normalize the data by the magnitude of the largest absolute anomaly magnitude in each window to avoid the possibility that any single eruption would dominate the epochal signal. Each window includes 11 years with 5 years before and 6 years after each eruption.

The average of the all 10 “full forcing” and 5 “volcanic only” CESM-LMEs is obtained first. For the multi-model ensemble (MME) mean, each model has the same composite weight, namely, the MME is calculated as the average of these 11 models after interpolating all model fields to 1° × 1° grid. The composite with each of the all “full forcing” and five “volcanic only” CESM-LMEs is performed for reference. We assess statistical significance of the MME at the 95% confidence level using Student’s t test, assuming that each model represents an independent sample for the MME.

3 Results

3.1 Observed and simulated volcano-El Niño relationship

In the instrumental observations for the 1963 Agung, 1982 El Chichón, and 1991 Pinatubo eruptions, an El Niño peaked in the first winter after each of the eruptions (Fig. 1). After these three tropical eruptions, WCEP westerly anomalies appeared before the peak phase of El Niño. Based on only three cases, we cannot conclude that these westerly anomalies are caused by tropical eruptions, which could be due to internal variability (Robock et al. 1995; Self et al. 1997).

Fig. 1
figure 1

Observed and simulated volcano-El Niño relationship. Five-month running average of Niño 3.4 index (K; dotted line: HadSST1; solid line: HadSST2; dot-dashed line: ERSSTv5) and 850 hPa WCEP (140° E–180°, 5° S–5° N) zonal wind (m s−1; NCEP-R1) anomalies in observations and reanalysis and in SEA composite of 41 ensemble members of CESM-LE (dashed line) for each of the 1963 Agung, 1982 El Chichón and 1991 Pinatubo eruptions. Solid segment of the curve denotes simulated anomaly significant above the 95% confidence level. Monthly zonal-mean AOD before and after the eruptions for these three eruptions are obtained from Ammann (2003). Year (0) denotes the eruption year, and year (1) is the first year after the eruption

Full size image

The SEA composite of the 41 CESM-LE ensemble members, after the internal mode is removed by using large ensemble mean, reproduces an El Niño response in the second winter after each of these three eruptions (Fig. 1). Before this peak phase of the El Niño, a significant WCEP westerly anomaly is already simulated in the eruption year for each eruption. Then, the El Niño develops and peaks in the second boreal winter. A significant positive Niño 3.4 index is simulated in the mean of these 41 CESM-LE ensemble members before the 1991 Pinatubo eruption, which means that even 41 ensemble members are not enough to totally remove a signal of internal variability. Results are not changed when using the original monthly data instead of smoothing data (Fig. S4). Similar results are also obtained at 1000 hPa (Fig. S5), except for a weaker amplitude than that at 850 hPa.

To identify the direct effect of volcanic forcing, the boreal fall is studied, after which the El Niño should mature in the boreal winter and the delayed responses are large. Before the mature phase of the simulated El Niño after 1963 Agung, 1982 El Chichón, and 1991 Pinatubo eruptions, positive SST anomaly is simulated in the Niño 3.4 region in the boreal fall (Fig. 2a) after each peak of these three eruptions (Fig. 1), although it is not significant. Associated with a significant WCEP westerly anomaly (Fig. 2a), significant dipole precipitation anomalies are also simulated (Fig. 2b), with a negative anomaly over the SSEA and a positive one over the WCEP. Significant negative precipitation anomalies are also simulated over the West African monsoon and ITCZ regions. Each of these three eruptions presents a similar significant dipole precipitation anomaly (Fig. S6). In the CESM-LE, the strength and meridional distribution of the 1962 Agung and 1982 El Chichón aerosols are similar. The stronger global cooling after 1962 Agung eruption than 1982 El Chichón eruption may be caused by global warming signal in 1982 response and different initial mean state. This is an interesting topic and will be studied in future. Our results doesn’t change when the global warming trend is removed (not shown). In the boreal summer, the WCEP westerly and dipole precipitation anomalies already occur, although the latter are not significant (Fig. S7). The first boreal winter has similar results as the boreal fall except that the positive SST anomaly in the Niño 3.4 region becomes significant (Fig. S8).

Fig. 2
figure 2

Simulated tropical anomalies after three recent large eruptions. a SST (shading; K) and 850 hPa wind (vector; m s−1) anomalies and b precipitation (shading; mm day−1) anomalies during September–November in the eruption year in SEA composite of 41 ensemble members of CESM-LE for the average of the 1963 Agung, 1982 El Chichón, and 1991 Pinatubo eruptions. Both stippling and vector denote anomalies significant above the 95% confidence level

Full size image

Based on the result that a significant WCEP westerly anomaly leads the mature phase of an El Niño, we conclude that there is a direct significant equatorial westerly response to a tropical eruption in the CESM-LE. The El Niño then develops in the second year through the Bjerknes feedback. The significant dipole precipitation anomalies over the SSEA and central equatorial Pacific may be one reason for this westerly anomaly excitation under large tropical eruption in the CESM-LE. Next, we look into this possibility by using more samples, namely, by using model outputs during the last millennium.

3.2 Westerly anomaly and precipitation suppression response

To check whether a westerly response exists in all models, we study the last millennium simulations by the PMIP3 and CESM-LME. In these simulations, the volcanic forcing in two different datasets usually peaks in the middle of the year, and the volcanic forcing from Gao et al. (2008) is usually larger than that from Crowley et al. (2008) (Fig. 3a). After the 10 large tropical eruptions, a significant decrease in tropical-mean net surface shortwave flux (Fig. 3b) and tropical-mean surface temperature (Fig. 3c) is simulated in all 10 PMIP3 models and in the multi-ensemble mean of the CESM-LME simulations. For the same volcanic forcing, the 10 PMIP3 models and CESM-LME have different radiative perturbations, and CCSM4 has the strongest volcanic-induced cooling, consistent with previous finding (Atwood et al. 2016). These different sensitivities to external forcing are caused by models’ unique implementation of volcanic forcing (Timmreck 2012; Toohey et al. 2014). Tropical cooling has a phase lag of several months after an eruption, due to slow oceanic responses (Iles et al. 2013). This tropical cooling at surface, along with global cooling, lasts for 2–3 years, consistent with previous finding (Sear et al. 1987; Mann et al. 1998; Crowley 2000; Robock 2000; Hegerl et al. 2003; Timmreck 2012; Liu et al. 2016). Since the volcanic forcing of Crowley et al. (2008) appears a few months earlier than that of Gao et al. (2008) (Fig. 3a), the model responses also reflect such time difference (Fig. 3b and c).

Fig. 3
figure 3

Simulated tropical responses to 10 largest tropical eruptions during last millennium. a Monthly zonal-mean AOD before and after the eruptions for the 10 largest eruptions during the last millennium from reconstructions of Gao et al. (2008) (blue line) and Crowley et al. (2008) (red line). Shading denotes one standard deviation among the 10 eruptions. b SEA composite of five-month running average in tropical (20° S–20° N) mean surface net incoming all-sky shortwave flux (SNSWF; positive: downward; dashed lines; W m−2) for the 10 eruptions in each of 10 PMIP3 models and CESM-LME. Warm and cool colors denote responses under the volcanic forcing of Crowley et al. (2008) and Gao et al. (2008), respectively. Thick solid segment of the curve denotes anomaly significant above the 95% confidence level. c and d are the same as b, except for the responses of tropical (20°S-20°N) mean surface temperature (TMT; K) and WCEP (140° E–180°, 5° S–5° N) westerly (m s−1) anomalies, respectively

Full size image

Before the first winter when an El Niño should peak, as in the observations and reconstructions (Khodri et al. 2017; Liu et al. 2018a, b), a significant WCEP westerly anomaly is excited in most of these 10 PMIP3 models and CESM-LME under a volcano-induced tropical cooling; the exceptions are the IPSL, BCC, and FGOALS models (Fig. 3d). The latter two models produce WCEP westerly anomalies, but they are not significant. A significant westerly anomaly in the MPI model is simulated before June since the strong forcing of Crowley et al. (2008) appears at the beginning of the eruption year (Fig. 3a). These significant WCEP westerly anomalies can even last till the second winter after the eruption in the CESM-LME. Since eight out of the 11 models simulate significant westerly anomalies, we conclude that there is a robust WCEP westerly response to tropical eruptions in these models.

The difference in AOD of Gao et al. (2008) and AOD of Crowley et al. (2008) is quite large (Fig. 3a), especially for the 1258 Samalas eruption (Table 2). Due to the effect of eruption strength on ENSO development (Emile-Geay et al. 2008; Lim et al. 2016), we perform a composite without the Samalas eruption, and the results do not change qualitatively (Fig. S9). These significant WCEP westerly anomalies can also be obtained when normalized SEA are performed to avoid the possibility that any single eruption would dominate the epochal signal (Fig. S10).

In a few months after a tropical eruption, a significant cooling is simulated mainly over all the continents along 20° S and 30° N, in the MME of the LME simulations (Fig. 4a). The cooling over the ocean is relatively weak; even small positive SST anomaly appears over the central equatorial Pacific. Associated with this zonal SST gradient, a significant WCEP westerly anomaly is simulated. Stronger cooling over the MC than its surrounding ocean, which is suggested to excite a westerly anomaly over the WCEP (Ohba et al. 2013), cannot be properly reflected in the MME of these models.

Fig. 4
figure 4

Precipitation suppression responses to 10 largest tropical eruptions during last millennium. a SEA composite of surface temperature (ST; shading; K) and 850 hPa wind (vector; m s−1) anomalies, b upward surface latent heat flux anomalies (shading; W m−2), and c precipitation anomalies (shading; mm day−1) and climatological precipitation (contour; mm day−1) for June–November averages of the 10 eruption years in the 10 PMIP3 models and in the CESM-LME. Multi-ensemble mean of the CESM-LME is performed before the composite; thus, each model has the same composite weight. Both shading and black vector denote anomalies statistically significant above the 95% confidence level. In c, climatological precipitation is contoured at 3 and 8 mm day−1

Full size image

Under the volcano-induced tropical cooling, the evaporation over all the tropics except the central-to-eastern equatorial Pacific is significantly reduced (Fig. 4b), which tends to reduce precipitation. After these 10 large tropical eruptions, energy budget analysis shows that the global Earth system loses heat due to reduction of irradiance, resulting in reduced global evaporation and precipitation, as well as less outgoing longwave radiation, and such hydrological weakening mainly occurs in the tropics (Fig. 5).

Fig. 5
figure 5

Global energy budget after 10 largest tropical eruptions during last millennium. Global and tropical (30° S–30° N) precipitation, latent heat flux and sum of TOA longwave, surface longwave, TOA shortwave, surface shortwave, and sensible heat fluxes (W m−2). Each term of the fluxes is also shown. Gray asterisks denote each model’s result. Positive value denotes upward at the TOA and downward at the surface

Full size image

Significant positive precipitation anomalies over the central equatorial Pacific and negative precipitation anomalies over the other tropical regions are seen in the simulations (Fig. 4c). Significant precipitation decrease only occurs over the climatologically wet regions, including the SSEA, ITCZ, and West African monsoon regions (Fig. 4c). The original, climatologically dry regions, especially the central-to-eastern equatorial Pacific, show no negative precipitation anomaly. Due to the strongest response occurs in the CCSM4 (Fig. 3), we also perform the composite without the CCSM4 model; similar results are obtained (Fig. S11). Results are not affected when removing the CESM-LME due to their large ensembles (Fig. S12), and the composite with each of the CESM-LME simulations also presents a similar result (Fig. S13).

Each of these 11 models presents a similar result (Fig. 6). Significant negative precipitation anomalies over the SSEA, ITCZ, and West African monsoon regions are simulated by all these models, although their strengths differ, with weak negative anomalies in BCC (Fig. 6e) and MIROC (Fig. 6k). The significant positive precipitation anomaly over the central equatorial Pacific, however, is only simulated by four (Fig. 6a, d, i, and j) out of the 11 models.

Fig. 6
figure 6

Precipitation suppression responses to 10 largest tropical eruptions during last millennium in individual models. SEA composite of precipitation anomalies (shading; mm day−1) and climatological precipitation (contour; mm day−1) for June-November averages of the 10 eruption years in each of the 10 PMIP3 models and in the CESM-LME. Multi-ensemble mean of the CESM-LME is performed before the composite; thus, each model has the same composite weight. Only precipitation anomalies significant above the 95% confidence level are shaded, and climatological precipitation is contoured at 3 and 8 mm day−1

Full size image

In summary, under tropical volcanic eruption, the tropics exhibits a precipitation suppression response over the SSEA, ITCZ, and West African monsoon regions due to reduction of evaporation. This is a “wet gets drier” response, which is first used to represent that the area experiencing drying generally fits very well with the climatologically wet region after large volcano eruptions (Iles et al. 2013).

3.3 Roles of different regional land cooling

Because of significant cooling over the extratropical continents (Fig. 4a), we further explore the role of different land cooling in the WCEP westerly anomaly excitation, and carry out seven sets of sensitivity experiments using the CESM1.2.2, and each set has 10 ensembles. The volcanic forcing-induced land cooling can be mimicked by increasing the land surface albedo. Associated with the global land cooling (Fig. 7a), suppressed precipitation is simulated over the SSEA, West African monsoon, and tropical American regions (Fig. 7b). The tropical land surface cooling is much weaker than extratropical land surface cooling due to cloud feedback, as the case done by the volcanic forcing (Fig. 6). Different from volcanic forcing-induced tropical cooling and precipitation suppression, the land-only cooling induces large positive precipitation anomalies over the tropical Indian Ocean, which could be a response to the westerly anomaly induced by the cooling of the West African monsoon regions. Significant WCEP westerly anomalies are simulated, along with low-pressure anomalies over the Pacific and Indian oceans and high-pressure anomalies over the land (Fig. 7c), consistent with previous finding for volcano-induced sea level pressure anomalies (Pausata et al. 2020).

Fig. 7
figure 7

Tropical circulation responses to surface cooling. Surface temperature (ST; shading; K), precipitation (shading; mm day−1) and sea level pressure (shading; hPa) with 850 hPa wind (vector; m s−1) anomalies during June-November responses to the surface cooling over the global (a–c) and extratropical (0–360°, 10°–90° S, 30°–90° N; d–f) regions. Both stippling and vector denote anomalies statistically significant at the 95% confidence level

Full size image

In the simulation with extratropical land cooling (Fig. 7d), there is no tropical precipitation suppression (Fig. 7e), and no significant WCEP westerly anomaly (Fig. 7f), although the mid-latitude land-sea thermal contrast and pressure gradient associated with extratropical continent cooling are present. Figure 8 shows the role of different regional cooling in the WCEP westerly anomaly excitation. The WCEP westerly anomaly is mainly excited by the tropical land cooling rather than by the extratropical continent surface cooling. The SSEA cooling and resulted precipitation suppression contribute more to this WCEP westerly anomaly than the tropical African cooling does. In the SSEA region, the WCEP westerly anomaly excited by the MC cooling is larger than that by AM cooling. Details of these responses can be found in Fig. S14.

Fig. 8
figure 8

Roles of different regions on WCEP westerly anomaly excitation. The WCEP (140° E–180°, 5° S–5° N) westerly anomaly caused by cooling over the global, extratropical (0–360°, 10°–90° S, 30°–90° N), tropical (0–360°, 10° S–30° N), SSEA (70–154° E, 10° S–30° N), AM (70–130° E, 10°–30° N), MC (96–154° E, 10° S–10° N), and tropical African (0–52° E, 10° S–30° N) regions. Black bars denote anomalies significant above the 95% confidence level

Full size image

These simulations confirm that the WCEP westerly anomaly is excited by tropical cooling-induced precipitation suppression, and SSEA cooling contributes to this anomaly more than the tropical Africa does. Both experiments with doubled albedo over the SSEA or tropical Africa show weak local land surface cooling and large precipitation suppression, similar as those shown in Fig. 7 for whole tropical land cooling, even though the tropics has stronger volcanic forcing than the extratropical regions. The tropical and extratropical responses to irradiance reduction are represented by precipitation suppression and land surface cooling, respectively (Figs. 4 and 7). Our results also demonstrate that for the MC cooling effect (Fig. 8) exciting the westerly anomaly (Ohba et al. 2013), it is the MC precipitation reduction, rather than land surface cooling (Figs. 4 and 7), that works.

3.4 Roles of different regional precipitation anomalies

To check whether a WCEP westerly anomaly can be excited by this precipitation suppression over the SSEA, ITCZ, and West African monsoon regions and to isolate the role of each of these precipitation anomalies, the Gill (1980) model is used. The Gill responses to the simulated precipitation anomaly after volcanic eruption can be obtained by specifying the condensation heating in Eq. (1). Following the relationship in Wang (1988), the non-dimensional condensational heating rate associated with precipitation is \(PL_{c} g/\Delta pC_{p}\), where \(L_{c} = 2.5 \times 10^{6}\, {\text{Jkg}}^{ - 1}\) is the latent heat of condensation, \(g = 9.8\, {\text{ms}}^{ - 2}\) is the gravity, and \(C_{p} = 1004 \;{\text{J kg}}^{ - 1}\, {\text{K}}^{ - 1}\) is the specific heat at constant pressure. The height scale of troposphere is \(\Delta p = 80,000{/}\pi \;{\text{Pa}}\) for the first baroclinic mode which has a vertical distribution of \({\sin}\left( z \right)\), where z changes from 0 to \({\uppi }\) according to 900 and 100 hPa, respectively. Precipitation anomalies \(P\) are specified by using the anomalies obtained from the average of the 11 models in Fig. 4c. The precipitation suppression response to tropical volcano is a common feature in all the models used here, while the positive precipitation anomaly over the central equatorial Pacific is not (Fig. 6). The positive anomaly over the central Pacific could also be a response to a WCEP westerly anomaly. Thus, in this work, we only use the negative tropical precipitation anomaly as the forcing in the Gill model.

Figure 9 shows the Gill responses to different volcano-induced tropical precipitation anomalies simulated by the 10 PMIP3 models and CESM-LME. The Gill responses can reproduce the circulation responses in these PMIP3 models and CESM-LME very well (Fig. 9a). Associated with these negative tropical precipitation anomalies, a strong WCEP westerly anomaly is excited in the Gill model, along with an easterly wind anomaly over the northern tropical Indian Ocean and westerly wind anomaly over the southern tropical Indian Ocean. The West African monsoon suppression can excite a strong Kelvin wave response to its east, and a westerly anomaly reaches the WCEP (Fig. 9b), which is consistent with previous finding (Khodri et al. 2017). A much stronger WCEP westerly anomaly is excited by the negative precipitation anomaly of the SSEA when compared to that by the West African monsoon (Fig. 9c). Similar result is obtained when the AM is removed since it excites a subtropical anticyclone and contributes to 28% of the westerly anomaly (Fig. S15), consistent with the sensitivity experiments (Fig. 8). The negative precipitation anomaly in the ITCZ, however, excites an easterly anomaly over the WCEP (Fig. 9d). The WCEP westerly anomaly is enhanced when the positive precipitation anomaly over the central equatorial Pacific, as shown in Fig. 4, is also included (Fig. S16).

Fig. 9
figure 9

Gill responses to volcanic forcing-induced precipitation anomalies over different regions. Shown are the lower-tropospheric wind (vector; m s−1) responses to the negative precipitation anomalies (shading; mm day−1) over a all tropics, b West African monsoon regions, c SSEA, and d ITCZ in Fig. 4c. A damping of 3.5 days is used in the Gill model

Full size image

Figure 10 shows the roles of different regional precipitation anomalies on a WCEP westerly anomaly excitation. The SSEA precipitation decrease accounts for 94% of a westerly anomaly excitation, while the West African monsoon suppression accounts for 49% of this westerly anomaly excitation. The negative precipitation anomaly in the ITCZ, however, tends to reduce 43% of this anomaly. The Gill response is sensitive to the damping scale. The West African monsoon suppression contributes more to this WCEP westerly anomaly excitation when a Kelvin wave is less damped with a damping of 15 days than of 3.5 days (Fig. 10).

Fig. 10
figure 10

Roles of different regions on WCEP westerly anomaly excitation. Ratio of WCEP (140° E–180°, 5° S–5° N) westerly anomaly caused by West African monsoon regions, SSEA, and ITCZ negative precipitation anomalies to that caused by full negative tropical precipitation anomalies. The blue and red bars denote results with a damping scale of 3.5 and 15 days, respectively, in the Gill model

Full size image

In conclusion, precipitation anomalies over the tropical Pacific are able to trigger the WCEP westerly anomaly. The negative precipitation anomalies over the West African monsoon regions and the SSEA all contribute to this westerly anomaly excitation, while the anomaly over the ITCZ does not. The SSEA contributes more to this westerly anomaly excitation than the West African monsoon regions do.

Here, we want to highlight the precipitation suppression over the SSEA, ITCZ, and West African monsoon regions since the mean-state-determined precipitation suppression without a negative precipitation anomaly over the central-to-eastern equatorial Pacific is important for the excitation of a WCEP westerly anomaly. We further perform an aqua-planet experiment in which the diabatic heating associated with the precipitation anomalies is zonally uniform. In this aqua-planet experiment, the zonally uniform tropical precipitation suppression after the tropical eruption cannot excite this type of westerly anomaly, and the weak tropical westerly anomaly is only excited by the quasi-geostrophical balance, as shown by the Gill model (Fig. 11).

Fig. 11
figure 11

Gill-type responses to zonally averaged precipitation anomalies. Shown are the lower-tropospheric wind (vector; m s−1) responses to a negative only tropical precipitation anomalies (shading; mm day−1) in Fig. 4c and to b zonally averaged precipitation anomalies. A damping of 3.5 days is used in the Gill model

Full size image

3.5 Failure to simulate El Niño in some models

Following the significant WCEP westerly anomaly an El Niño should develop via the Bjerknes feedback (Bjerknes 1969). The relative SST anomaly is used to identify an El Niño response under overall tropical cooling (Maher et al. 2015; Khodri et al. 2017), which is obtained by removing the tropical-mean (20° S–20° N) SST anomaly. The simulated relative Niño 3.4 in the first boreal winter after an eruption increases with the preceding WCEP westerly anomaly, with a linear correlation of 0.61 (p < 0.05) among these different models (Fig. 12a). However, three models, i.e., BCC, MPI, and MIROC, fail to simulate a positive relative Niño 3.4 index after a WCEP westerly anomaly. The underestimated Bjerknes feedback, found in both the BCC and MPI (Bellenger et al. 2014), is extreme and may limit their ability to reproduce the positive relative Niño 3.4 index.

Fig. 12
figure 12

Failure to simulate El Niño in some models. a Scatter plots of SEA composite of December–February relative Niño 3.4 index anomalies vs preceding June-November WCEP westerly anomalies during the eruption year for each of the 10 PMIP3 models and for the CESM-LME. Multi-ensemble mean of the CESM-LME is performed for comparison. Asterisks denote significant at the 95% confidence level. b and c are the same as a, except for the relative Niño 3.4 index anomaly vs central equatorial Pacific (160° E–160° W, 5° S–5° N) precipitation anomaly and for Niño 3.4 index vs WCEP westerly anomaly, respectively

Full size image

Although most of the models can simulate the precipitation suppression response and associated WCEP westerly anomaly, only a few of them can simulate the positive precipitation anomaly over the central equatorial Pacific (Fig. 6). The simulated relative Niño 3.4 index also increases with the preceding positive precipitation anomaly over the central equatorial Pacific, with a linear correlation of 0.73 (p < 0.05) among these different models (Fig. 12b). This result means that the models with more sensitive convective feedback can simulate an El Niño after the eruption more easily. This result is also consistent with previous finding for ENSO simulation (Neale et al. 2008; Guilyardi et al. 2009; Watanabe et al. 2011a). Among the linear atmospheric feedbacks driving ENSO evolution (Jin et al. 2006), the Bjerknes feedback (Bjerknes 1969) is still strongly underestimated by most models, which is sensitive to the convective scheme (Watanabe et al. 2011a).

In both observations (Khodri et al. 2017) and reconstructions (Liu et al. 2018a, b), a positive Niño 3.4 index is found in the first winter after tropical eruptions. Only two models (FGOALS and MRI) can simulate positive Niño 3.4 index, and nine models fail to simulate El Niño after a WCEP westerly anomaly (Fig. 12c). The El Niño can be simulated in the second winter by the CESM-LME. In the CCSM4, the central-to-eastern Pacific is dominated by strong volcanic cooling, which suppresses the Bjerknes feedback-induced warming. This difference between SST and relative SST anomalies is also found in other models (Khodri et al. 2017), and the failed simulation of the El Niño in the first boreal winter after the eruption is in contrast to the observation (Khodri et al. 2017) and reconstructions (Liu et al. 2018a,b).

In summary, a WCEP westerly anomaly alone does not seem to be sufficient for a realistic simulation of further development of El Niño, and the failure to simulate El Niño can also occur when the WCEP westerly anomaly is simulated. The models with more sensitive convective feedback can simulate El Niño more easily.

4 Concluding remarks

In this study, we explore the mechanism for a westerly anomaly excitation over the WCEP in multiple models. The CESM westerly anomaly should favor the El Niño development. We conclude that there is a robust WCEP westerly response to the large tropical eruption in these models, based on the finding that eight out of the 11 studied models can simulate a significant WCEP westerly anomaly after a large tropical eruption. Within a few months after the large tropical eruption, tropical evaporation is decreased due to a reduction of incoming shortwave radiation by stratospheric sulfate aerosols. Under this overall tropical drying, there are negative precipitation anomalies over SSEA, West African monsoon, and ITCZ regions. Sensitivity experiments by the CESM show that the WCEP westerly anomaly can be excited by the tropical land cooling-induced precipitation suppression, rather than by the extratropical land surface cooling, and SSEA cooling contributes more than the West African monsoon does in terms of exciting the WCEP westerly anomaly. Theoretical results show that a WCEP westerly anomaly is excited due to a Gill response to weakened precipitation over the SSEA and West African monsoon regions, with the SSEA contributing more to this westerly anomaly excitation than the West African monsoon does, while the ITCZ weakening excites an easterly wind anomaly over the WCEP.

The simulation of WCEP westerly anomaly does not guarantee the simulation of El Niño in some models, and the simulation of El Niño can also fail. The models with more sensitive convective feedback tend to simulate an El Niño more easily. The upcoming Model Intercomparison Project on the climate response to volcanic forcing (VolMIP) (Zanchettin et al. 2016), with more models participating, will allow us to validate the results obtained here using more models. A priority for future studies should be investigating the failure to simulate the El Niño in some models, after a robust westerly anomaly is reproduced. This calls for further studies on model sensitivity to radiative forcing.