Statistical testing in observations

We investigated the equatorial Pacific (Niño3.4 and Niño3 regions) SST response to major tropical volcanic eruptions of the instrumental period using composites in HadISST observations46. We restricted our analyses to 1870–2010 as there are too many gaps in spatial sampling before 1870. The earliest part of the record (late nineteenth century) is, however, less reliable due to lower spatial resolution. We build composites using the five tropical volcanic eruptions in 1870–2010 with significant injection of sulphuric aerosol in the tropical stratosphere: Pinatubo (June 1991), El Chichón (April 1982), Mt Agung (March 1963), Santa María (October 1902) and Krakatau (August 1883). SST anomalies are computed relative to the 5-year mean preceding each eruption (to remove multidecadal variability and long-term trends) and smoothed with a 3-month Hanning filter (to filter out intraseasonal variability). To evaluate the level of statistical significance, we first count the number of events with anomalies that have the same sign, as shown in Fig. 1a. If we assume that there is the same probability of having positive and negative signs, according to the binomial law, the likelihood for having 5 anomalies of the same sign is 0.03, i.e., corresponding to a 97% confidence level. Second, we use a Monte Carlo approach to compute the distribution for the composite value of November–January RSST average anomalies of 5 randomly picked years within the historical period, the result of which is given in Supplementary Fig. 2. The probability in the Niño3.4 composite of RSST anomaly equal or above the observed 1.1 °C value is very small (0.3%). The level of statistical significance is similar in other data sets such as HadiSSTv1.1, HadSST31, ERSSTv2b247, Kaplan SST348 and Niño34hist449. Various methods to compute the anomalies (15-year high-pass filter, removing the mean value for the previous 1, 3, 5 or 10 years) were tested and did not affect the results.

CMIP5 historical simulations

We considered 26 CMIP5 models and a total of 106 members from the CMIP5 historical simulation database31 (Supplementary Table 1). Each member is a coupled ocean–atmosphere simulation using the CMIP5 recommended transient forcing from 1850–2005. We did not include simulations using anomalies of the solar constant to mimic volcanic forcing, interactive chemistry, carbon cycle or prognostic global aerosol models in our analyses, to restrict the potential noise associated with those extra degrees of freedom. Despite this, there is still a relatively large diversity of specified external forcings, as shown in Supplementary Table 1. Anomalies were computed relative to the 1961–1990 climatology of each model simulation. We analysed SST evolution in these simulations by using both the raw and relative (obtained by removing the global tropical 20 °N–20 °S mean SST) SST anomalies. For each of the 106 CMIP5 historical simulations31, we identified El Niño events from averaged Niño3.4 region monthly raw and relative SST standardized anomalies that exceed half a standard deviation over the 1870–2010 simulated period. The percentage of CMIP5 simulations displaying an El Niño-like state and its statistical significance are shown in Supplementary Figs. 3–5. Using various methods to compute the SST anomalies (removing the mean climatology of the previous 5 years, using the 1961–1990 or the 1870–2010 baseline climatology) do not affect the results. The difference between raw SST (Supplementary Fig. 1a, b) and RSST (Fig. 1a, b) is weak in observations, but it is large in CMIP5 models (see Supplementary Fig. 5). The largest uncertainties in the estimates of radiative forcing from CMIP5 historical simulations occur during periods of volcanic activity43 and generally models tend to overestimate the observed post-eruption global surface cooling44. Relying on relative SST and precipitation anomalies partly corrects this bias, so that it reveals the dynamical ENSO response to volcanic forcing in CMIP5 models.

Coupled model simulations for the 1991 Pinatubo eruption

To help analyse the response of ENSO to volcanic forcing, a suite of short ensemble experiments was run over January 1991 to December 1993 for the Mt Pinatubo eruption (June 1991) using the IPSL-CM5B-LR28 ocean atmosphere CGCM. The stratospheric aerosol optical depth (hereafter AOD) forcing due to the Pinatubo eruption in each member uses the Gao et al.51 data set. The stratospheric volcanic cloud spreads from the tropics to the poles from the start of the eruption in June 1991 until late 1992, and then slowly disappears as the aerosols are totally washed out of the atmosphere in 1993. The total AOD reaches its maximum in boreal autumn–winter 1991–1992. Corresponding control ensembles starting from the same initial conditions but without Pinatubo volcanic forcing were performed.

IPSL-CM5B-LR coupled model simulations ensemble

We use the IPSL-CM5B-LR CGCM. The Laboratoire de Météorologie Dynamique (LMD) developed the LMDZ5B atmospheric component, which is described in Hourdin et al.52. This model uses hybrid σ–p vertical coordinates, with 39 pressure levels including 18 in the stratosphere. The LMDZ5B model uses ORCHIDEE as land surface component. We used the low-resolution grid with 1.87° in latitude and 3.75° longitude. The oceanic component, Nucleus for European Modeling of the Ocean (NEMO) version 3.2, is an oceanic GCM with 31 vertical levels (whose thickness varies from 10 m near the surface to 500 m towards the bottom), a 2-level representation of sea ice and a mean spatial horizontal resolution of about 2° (with a refinement of the latitudinal resolution to 0.5° near the equator; ORCA2 grid). Previous works53 found that this model reproduces ENSO variability with two spectral peaks around 3–3.5 years and beyond 4 years, which is in good qualitative agreement with observations.

The observed ENSO seasonal phase locking with a peak occurring mostly in boreal winter (November to January) and feedbacks are also well captured in this model53. Initial conditions were chosen in the IPSL-CM5B-LR standard historical run. To disentangle the role of Pinatubo eruption on ENSO evolution and avoid model spin-up, we chose restart days in the historical run with greenhouses gases and tropospheric aerosols concentrations close to the 1990s levels and more than 5 years after any volcanic eruption. We selected three restart dates corresponding to 1 June (i.e., shortly before the Pinatubo eruption) with initial conditions that would be favourable to the development of an El Niño, neutral, or La Niña phase in the absence of volcanic forcing. WWV (the volume of water above the 20 °C isotherm within 5 °N–5 °N 120 °E–80 °W) is a good precursor of the upcoming ENSO phase in nature4. This is also the case in our model as demonstrated by a peak correlation of r = 0.7 when the Niño 3.4 SST index lags the WWV in May by 6 months in the IPSL-CM5B CMIP5 historical simulation54. We hence selected 3 model states representing discharged, near-neutral and recharged WWV based on normalized WWV values in May. For the three initial states, 30 members (forced and unforced) were generated by adding a small white noise on the initial SST. Consistent with observations, the discharged, neutral and recharged states ensemble mean respectively produce a La Niña, near-neutral and El Niño state in the absence of volcanic forcing (see Fig. 2 and Supplementary Fig. 7).

CNRM-CM5 coupled model simulations ensemble

To evaluate the robustness of the results obtained with the IPSL climate model, we also used the CNRM-CM5 model developed jointly by CNRM-GAME (Centre National de Recherches Météorologiques—Groupe d’études de l’Atmosphère Météorologique) and Cerfacs (Centre Européen de Recherche et de Formation Avancée). CNRM-CM5, described by Voldoire et al.50, couples the ARPEGE-Climat (v5.2) atmosphere, NEMO (v3.2) ocean, the ISBA land surface and GELATO (v5) sea ice models through the OASIS (v3) coupler. The horizontal resolution is 1.4° in the atmosphere and 1° in the ocean (with a refinement of the latitudinal resolution to 0.5° near the equator). The model reproduces ENSO typical variability well with an ∼2‒7-year period and a seasonal phase locking with a peak occurring mostly in boreal winter, in good qualitative agreement with observations53. The ensemble of the CNRM-CM5 model was designed to explore the role of Pinatubo eruption on ENSO response when starting from random initial conditions. The initial dates were chosen randomly in the control and historical runs performed with the same model version50.

The two-tier sensitivity experiments approach

To explore the role of various processes that drive the climate response to volcanism, we ran two-tier experiments using an atmospheric model simulations ensemble (AGCM) and a linear ocean model.

The atmospheric model simulations ensemble

We ran several 30-member AGCM ensembles using the atmospheric component of the IPSL-CM5B coupled model52. Daily fields from each of the 30 members of the coupled model ensembles starting from neutral initial conditions are used to derive a set of boundary conditions (described below) for the stand-alone atmospheric component (LMDZ AGCM) of the IPSL-CM5B model, which includes the ORCHIDEE land surface module. The stand-alone atmospheric model is run for 2 years from 1 January 1991 onward (i.e., year 0; 5.5 months before the Pinatubo eruption) to enable the atmosphere to dynamically adjust to the boundary conditions. We have first run a Control 30-member AGCM ensemble using surface boundary conditions (i.e., SST and sea ice cover) from the coupled model control ensemble without including the changes induced the volcanic aerosol. The ALL experiment follows the same protocol, but includes the stratospheric AOD evolution corresponding to the Pinatubo eruption as in the coupled model experiments and specified SST and sea ice cover from each member of the coupled model Pinatubo ensemble. The anomalies of ALL relative to the Control are almost identical when compared to the anomalies of the coupled model experiment (Supplementary Fig. 10). The small differences between the coupled and the atmospheric experiment are explained by the differences in the surface flux due to coupling55. To address the role of the atmospheric, land and ocean responses in triggering surface wind anomalies, five complementary sets of experiments have been run (see summary in Supplementary Table 2).

The ATM, LAND and OCEAN experiments were respectively designed to represent the effect of aerosol forcing through its direct impact on the atmospheric structure (ATM), and indirect impacts through changes in land (LAND) or ocean (OCEAN) surface temperatures. The ATM and LAND scenarios were developed using the same protocol as in ALL but with prescribed daily SST from the Control members (i.e., with a SST that excludes the impact of volcanic aerosol forcing). Volcanic aerosol forcing is included in the ATM experiments, but the surface albedo in this experiment is modified so that continental surfaces do not cool in response to volcanic forcing (details below). In LAND, there is no prescribed aerosol forcing, but the albedo modification enforces a land surface cooling that is consistent with that of the Pinatubo coupled model experiment. This strategy has the advantage of allowing surface land temperature to vary in the presence of SST anomalies (in the OCEAN experiment), as is seen in other model experiments in the context of global warming scenario37, 38. This is hence more physically relevant than experiments that constrain land surface temperature19 and which can result in unphysical and exaggerated land cooling, and hence exaggerated land–ocean gradients. We modify the land surface albedo to restore the atmospheric temperatures over land to that of the Control coupled model ensemble in ATM and to that of the Pinatubo coupled model ensemble in LAND.

The boundary layer temperature changes ΔT in response to a given change in the net shortwave radiation ΔSWnet SFC at the surface can be approximated assuming linear dependency such as:

$$\lambda \Delta T = \Delta \left( {SWdow{n_{{\rm{SFC}}}} - SWu{p_{{\rm{SFC}}}}} \right) = \Delta SWne{t_{{\rm{SFC}}}}$$

with λ being the feedback parameter including that due to the Planck emission (i.e., the feedback associated with changes in surface temperature), changes in vertical temperature lapse rate, water vapour and clouds. Here, we specify SWnet SFC by prescribing the surface albedo α. The surface albedo is related to the net shortwave radiation by:

$$SWne{t_{{\rm{SFC}}}} = SWdow{n_{{\rm{SFC}}}}\left( {1 - \alpha } \right)$$ (1)

In particular, we have:

$$SWne{t_{{\rm{cSFC}}}} = SWdow{n_{{\rm{cSFC}}}}\left( {1 - {\alpha _{\rm{c}}}} \right)$$ (2)

$$SWne{t_{{\rm{pSFC}}}} = SWdow{n_{{\rm{pSFC}}}}\left( {1 - {\alpha _{\rm{p}}}} \right)$$ (3)

Here, the subscript c denotes the values of the control coupled model experiment, while the subscript p designates the values of the Pinatubo coupled model experiments. Then, we calculate Δα, the relative change of the incoming SW radiation due to stratospheric sulphate aerosols:

$$\Delta \alpha = \left( {SWdow{n_{{\rm{pSFC}}}}-SWdow{n_{{\rm{cSFC}}}}} \right)/SWdow{n_{{\rm{cSFC}}}}$$ (4)

In ATM, the downward shortwave radiation is modified due to the direct radiative effect of stratospheric sulphate aerosols (i.e., SWdown SFCATM = SWdown pSFC ), but we modify the surface albedo to α ATM so that the surface-absorbed net shortwave radiation remains equal to that of the control run (i.e., SWnet SFCATM =SWnet cSFC ). Using previous relationships together with (2) and (4), we obtain:

$${\alpha _{{\rm{ATM}}}} = \left( {{\alpha _{\rm{c}}} + \Delta \alpha } \right)/\left( {1 + \Delta \alpha } \right)$$ (5)

Conversely, in LAND, the surface albedo is set to the α LAND value, with downward shortwave radiation equal to the control experiments (SWdown SFCLAND = SWdown cSFC ), but surface-absorbed shortwave equal to that of the Pinatubo experiments (i.e., SWnet SFCLAND =SWnet pSFC ):

$${\alpha _{{\rm{LAND}}}} = {\alpha _{\rm{p}}} - \Delta \alpha + {\alpha _{\rm{p}}}\Delta \alpha $$ (6)

We diagnosed Δα from the control and Pinatubo mean ensemble simulations and performed ATM and LAND simulations imposing the albedo in order to restore the surface temperature towards control and Pinatubo values respectively. This protocol involves a linearity assumption of assumed constant values for the feedback parameter λ. Surface temperature anomalies in the ATM, LAND and OCEAN experiments (Fig. 4 and Supplementary Fig. 10) however add up quite linearly to surface temperature anomalies in the ALL atmospheric or Pinatubo coupled model experiment, indicating the validity of our approach. The LAND experiment isolates the effect of volcanically induced continental cooling but does not narrow down the region that plays a key role for inducing wind anomalies over the tropical Pacific. We hence also developed a suite of additional 30-member LAND sensitivity experiments (with and without Pinatubo radiative forcing) that allows exploring the respective roles of land cooling in the tropics (LAND-T), the extra-tropics (LAND-ET), Africa (LAND-Africa), Southeast Asia (LAND-SEA) and the maritime continent (LAND-MC) only (See Supplementary Table 2).

The linear Indo-Pacific Ocean model

The atmospheric experiments above isolate the wind changes over the Pacific associated with three different processes: direct effect of radiative forcing on the atmosphere (ATM) and indirect effect through the cooling of continents (LAND) and induced oceanic SST gradients (OCEAN). To translate those wind signals into SST responses, we use a linear continuously stratified56 (LCS) ocean model of the Indo-Pacific region combined with a linear SST equation. The LCS model is forced by the atmospheric forcing produced by each of the 30-member experiments above to deduce the equatorial Pacific Ocean response (and in particular the SST response). The LCS model for the Indo-Pacific Ocean resolves vertical baroclinic modes. It is used to simulate thermocline depth and near-surface current anomalies. It is combined with a linear SST equation57 to realistically resolve the variations of SST due to wind stress anomalies:

$$\frac{{{\rm d}SST\prime }}{{{\rm d}t}} = - U\prime \frac{{{\rm d}\overline {SST} }}{{{\rm d}x}} - V\prime \frac{{{\rm d}\overline {SST} }}{{{\rm d}y}} + \gamma H\prime - \beta SST\prime $$

The horizontal advection term is approximated from the anomalous advection of the climatological SST (\(\overline {SST} \)) gradient (here from the IPSL historical run climatology) by the LCS 0–100 m averaged anomalous currents U' and V'. In the Indo-Pacific region, zonal advection dominates over meridional advection (not shown, see for example Vialard et al.35). The γH' term (H' is the thermocline depth anomaly obtained as H'=150 SSH' from the sea surface height anomaly SSH') parameterizes the vertical physics, i.e., the interannual modulation of mixing and upwelling by thermocline depth anomalies. The coefficient γ is constant to the east of 140 °W (γ=0.028 °C/month, as in Burgers58), and linearly decreases until 1/5 of its maximum value at the western boundary of the Pacific, as in McGregor et al.59. This mimics the decreasing effect of the upwelling as the climatological thermocline deepens towards the west. Heat fluxes are parameterized as a Newtonian damping –β SST′, with the associated damping timescale β−1=2 months, close to values used by McGregor et al.59 and Burgers58. The 5 first baroclinic modes are used (even though mainly only the first two modes contribute)57. This model has the advantage of being relatively realistic (Validation: high correlations of 0.8–0.9 with observations (TAO, AVISO, OSCAR), notably for zonal current and SST in Niño4 (central Pacific), for which the first two baroclinic modes are important, and also for SST and SSH in Niño3 (eastern Pacific) when the model is forced by ERAI (see details in Izumo et al.57), as well as rather simple and cheap to use. It is used to simulate the SST response to wind stress anomalies diagnosed from the ATM, LAND and OCEAN AGCM ensemble members. The high correlation of the IPSL-CM5B coupled model Niño3.4 SST anomalies with those simulated by the LCS model forced by the coupled model experiments (r = 0.9; Supplementary Fig. 9) and ALL (r = 0.8, not shown) ensemble wind stress justifies our approach. Furthermore, the favourable comparison of the ALL ensemble mean with the coupled model experiment wind stress anomalies in the western equatorial Pacific region (Supplementary Fig. 10) lends credence to our method and suggests that our two-tier methodology to investigate various physical mechanisms in forced experiments derived from the coupled one is sound.

Estimating signal-to-noise ratios

For the analysis of these Pinatubo sensitivity experiments, we used two approaches to evaluate the statistical significance of anomalies associated with volcanic forcing, relative to purely random, internal coupled ocean–atmosphere variability. The anomalies associated with the response to volcanism are computed as paired differences between each of the members of the Pinatubo experiment and its unforced counterpart. The first statistical calculation estimates the statistical significance of the anomaly between the two ensemble means using a two-tailed Welch’s t-test, and taking into account the size and variance of each (control and Pinatubo) ensemble. For CMIP5 model analyses, the significance of anomalies relative to the climatology is measured using a two-tailed Student’s t-test, considering each model as an independent degree of freedom. For both our paired Pinatubo ensembles and CMIP5 models, we also estimated at each time step and grid point the percentage for which at least two thirds of the individual simulations have anomalies of consistent sign, which is a commonly used technique for analysing outputs of CMIP models. This allows us to verify that statistically significant anomalies in the ensemble mean are representative of the ensemble members and not due to a few extreme members. Figure captions specify the test for each figure.

Code and data availability

The code and data that support the findings of this study are available from the corresponding authors on request.