We analyzed global surface air temperature (SAT) responses to five major tropical volcanic eruptions from 1870 to 2005 using outputs from 97 historical and 58 Atmospheric Model Intercomparison Project (AMIP) runs that participated in phase 5 of the Coupled Model Intercomparison Project (CMIP5). In observations, there was a 3-yr global cooling trend after the eruption due to reduced shortwave radiation, and a 0.1-K average global-mean SAT recovery, consisting of El Niño–like tropical warming and Eurasian warming, occurred in the first posteruption boreal winter. This global cooling pause was simulated by the multimodel ensemble (MME) mean of the AMIP runs, but not the MME of the historical runs due to the absence of El Niño–like warming. In the historical runs, simulation of El Niño–like warming was influenced by the initial ocean condition (IOC). An El Niño–like response was simulated when the IOC was not in an El Niño state, but the warming was much weaker compared to observations. The Eurasian warming response, despite being reproduced by the MME mean of both AMIP and historical runs, was not as strong as in observations. This is because the simulated positive polar vortex response, an important stratospheric forcing for Eurasian warming, was very weak, which suggests that the CMIP5 models, and even the Climate Forecast System model, underestimate volcanic effects on the stratosphere. Most of the coupled models failed to replicate both the El Niño and the enhanced polar vortex responses, indicating an urgent need for improving air–sea interaction and stratospheric processes in these models.
Large volcanic eruptions provide an important external forcing that has caused significant global surface cooling typically lasting for 3 years (Crowley 2000; Hegerl et al. 2003; Mann et al. 1998; Robock 2000; Sear et al. 1987; Timmreck 2012). This cooling process starts when a large amount of SO2 gas injected into the stratosphere reacts with OH and H2O to form sulfate aerosols, which can perturb the radiative balance (Robock 2000). For example, the year after the strong eruption of Tambora in 1815 is known as the year without a summer (Oppenheimer 2003). Surface cooling, however, is not always the case after a large tropical eruption. The volcanic eruption–induced global cooling was absent due to El Niño in the first posteruption boreal winter after the 1963 Agung eruption (Lehner et al. 2016) and the 1982 El Chichón eruption (Lehner et al. 2016; Santer et al. 2014), thereby forming a pause in the 3-yr global cooling following those tropical eruptions.
Changes in global mean surface air temperature (SAT) are controlled not only by external forcing but also by El Niño–Southern Oscillation (ENSO). Over the past decades, a positive global mean SAT anomaly is usually related to El Niño (Trenberth 2002; Trenberth and Fasullo 2013), which is the warm phase of ENSO. Global cooling emerges after a large volcanic eruption when the El Niño signal is removed from observations (Angell 1988; Gu and Adler 2011; Nicholls 1988; Thompson et al. 2009), from simulations (Paik and Min 2016), or from both observations and model simulations (Santer et al. 2014; Soden et al. 2002), which suggests that ENSO may be one reason for the inconsistence in global cooling in the first boreal winter after a large volcanic eruption.
El Niño occurred after four of the five large volcanic eruptions during 1870–2005 (Khodri et al. 2017). The relationship between El Niño events and volcanic eruptions in this period was originally considered a coincidence (Robock 2000; Self et al. 1997). However, analyses of long-term reconstruction data have suggested that large tropical eruptions can increase the likelihood of El Niño (Adams et al. 2003; Liu et al. 2018b), mainly through the dynamic thermostat mechanism (Clement et al. 1996; Maher et al. 2015; Mann et al. 2005; Ohba et al. 2013; Predybaylo et al. 2017) and also via other mechanisms as reviewed in Liu et al. (2018a).
Surface warming over the Eurasian continent was also observed to reach a maximum of 3 K in the first boreal winter after the 1991 Pinatubo eruption (Robock 2002; Robock and Mao 1992), which might also contribute to the global cooling pause following the eruption. The stratospheric pathway is a potential explanation for this boreal-winter warming after a large eruption. Volcanic aerosols from a large tropical eruption warm the low-latitude lower stratosphere; the resulting temperature gradient toward the equator decreases the high-latitude geopotential height and strengthens the stratospheric polar vortex (Robock 2000; Stenchikov et al. 2002). Positive westerly anomalies associated with enhanced polar vortex will trap tropospheric wave energy via planetary wave reflection (Graf et al. 2007; Perlwitz and Harnik 2003). In addition, the positive Arctic Oscillation or the North Atlantic Oscillation (NAO) dominates the boreal winter circulation, resulting in Eurasian winter warming from January to February (Graf et al. 2014; Perlwitz and Graf 1995).
It is not clear whether the state-of-the-art models in phase 5 of the Coupled Model Intercomparison Project (CMIP5) can simulate the pause of global cooling, including both tropical El Niño–like warming and Eurasian warming, during the first boreal winter of the 3-yr global cooling after a large tropical eruption. Although the models successfully reproduce El Niño responses in the Pacific sea surface height, their multimodel ensemble (MME) mean usually fails to simulate the observed response in sea surface temperature (SST) to the large tropical volcanic eruptions during the period of 1870–2005 (Ding et al. 2014; Maher et al. 2015). Focusing on the largest eruptions of Krakatau and Pinatubo, Zambri and Robock (2016) reported that some CMIP5 models replicate the observed warming over northern Europe and Asia and the observed polar vortex enhancement in the first posteruption boreal winter. When some relatively weaker eruptions are included, the CMIP5 models generally fail to capture the dynamic responses in the Northern Hemisphere (NH) after most tropical eruptions, such as the NAO, polar vortex, and Eurasian warming (Bittner et al. 2016; Driscoll et al. 2012).
In this work, we investigate the fidelity of SAT responses in CMIP5 models during the first boreal winter after a large tropical volcanic eruption. The paper is organized as follows. In section 2, we present the data and methods. Simulations of global cooling pause and El Niño response are discussed in section 3. The Eurasian SAT and polar vortex responses are the foci of section 4. In section 5, we discuss the internal variability–induced global cooling pause in the Community Earth System Model Large Ensemble (CESM-LE). Conclusions are given in section 6.
2. Data and methods
Monthly reanalysis data of geopotential height, air temperature, and zonal wind were obtained from 1) the Twentieth Century Reanalysis version 2c (20CRv2c; Compo et al. 2011) for the period of 1851–2014 and 2) the Japanese 55-year Reanalysis (JRA-55; Kobayashi et al. 2015) for the period of 1958–2013. The long-term observed monthly SAT field from 1870 to 2018 was also used; the data were from the 1200-km smoothed NASA Goddard Institute for Space Studies (GISS) Surface Temperature Analysis (GISTEMP) version 3 (Lenssen et al. 2019). Monthly SST data were from the following sources: 1) the Extended Reconstructed Sea Surface Temperature dataset version 5 (ERSSTv5; Huang et al. 2017) for the period of 1854–2017, 2) the Hadley Centre Sea Ice and Sea Surface Temperature dataset version 1 (HadISST1; Rayner et al. 2003) for the period of 1870–2016, and 3) the uninterpolated Second Hadley Centre Sea Surface Temperature dataset (HadSST2; Rayner et al. 2006) for the period of 1850–2004.
Five large tropical volcanic eruptions during 1870–2005 were selected in this study, namely, the eruption of Krakatau in August 1883, of Santa Maria in October 1902, of Agung in March 1963, of El Chichón in April 1982, and of Pinatubo in June 1991 (Table 1). They have the five largest aerosol optical depths (AOD) at 550 nm among the tropical volcanic eruptions since 1870 (Khodri et al. 2017; Sato et al. 1993). 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), was not considered in this work. Since the global cooling pause as well as the El Niño and Eurasian SAT responses peaked in the first posteruption boreal winter, we use calendar months rather than the months after the eruption peak in our study.
Table 2 gives details of the 32 coupled general climate models (GCMs) for historical runs covering the period of 1870–2005 (Taylor et al. 2012). The 17 atmospheric GCMs for the Atmospheric Model Intercomparison Project (AMIP) runs covering the period of 1979–2005, along with observed SSTs and sea ice (Taylor et al. 2012) used in the AMIP runs, are also introduced in the table. The models in the historical and AMIP runs were forced by both anthropogenic and natural forcing, including the five large tropical eruptions in the historical runs and the latest 1982 El Chichón and 1991 Pinatubo eruptions in the AMIP runs. Due to the important role of large ensemble members in detecting responses of ENSO (Maher et al. 2018) and NH polar vortex (Bittner et al. 2016) to volcanic forcing, 97 available ensemble members of the historical runs for the 32 models and 58 available ensemble members of the AMIP runs for the 17 AMIP models were considered. Table 2 also lists the volcanic forcing datasets used for these five eruptions, which were obtained from Sato et al. (1993), Andres and Kasgnoc (1998), Stenchikov et al. (1998), Ammann et al. (2003), and Ammann et al. (2007). Models adopting insolation reduction rather than aerosol forcing were not considered in this work. All model outputs were interpolated onto 2.5° × 2.5° grid using bilinear interpolation. For surface temperature, 2-m SAT was used.
In addition, we employed the fully coupled CESM-LE, which has 41 ensemble members at roughly 1° resolution for 1920–2005 following the historical run’s setup and contains three selected eruptions, namely the 1963 Agung, 1982 El Chichón, and 1991 Pinatubo eruptions. A full description of the CESM-LE can be found in Kay et al. (2015).
To isolate climate responses to the five eruptions from the background noise, we obtained anomalies by removing the average of climatology in the five years preceding each eruption. The long-term linear trend was also removed. Since the AMIP runs started from 1979, we used the average of climatology in the three years preceding the 1982 El Chichón eruption. When calculating the MME mean, averages were taken among ensemble members of each model first before getting the multimodel mean; in this way, all models were given equal weight. Similar results were also obtained by directly performing the average across ensemble members of all models without considering model bias (not shown).
Significance was calculated using the bootstrapped resampling method with 10 000 times of random synthesis across the whole period, and the 95% confidence level came from the 2.5%–97.5% range of the 10 000 draws. For the MME mean, we calculated the number of models that agreed on the sign of the response. In this way, the overall response could be considered without disturbances from large responses that may be dominated by only a few models (Barnes et al. 2016; Zambri and Robock 2016). If the data were purely random, for the historical runs the probability of at least 22 of the 32 models agreeing on the sign of the response was 5.0% with calculation based on the binomial distribution, which equaled the typical 95% confidence limit. For the AMIP runs, the probability of at least 13 of the 17 models agreeing on the sign of the response was 4.9%, which was close to the typical 95% confidence limit.
The Niño-3.4 index was defined as the SST anomaly averaged over the region of 5°S–5°N, 170°–120°W. We also used the relative SST anomaly reflecting the tendency of ENSO isolated from volcanically induced surface cooling (Khodri et al. 2017), which was obtained by removing the average of the tropical (20°S–20°N) SST anomaly from the original SST anomaly. To describe the stratospheric response in the North Pole region, the polar vortex index was defined as the inverted geopotential height anomaly averaged over the North Pole (65°–90°N, 0°–360°) at 50 hPa during the boreal winter. Here, boreal winter refers to December–February (DJF). Year (0) is defined as the eruption year, and year (−1) represents the year before the eruption.
3. Global cooling pause and El Niño response
Figure 1 shows the global mean SAT anomalies for the five large tropical volcanic eruptions and their composite values from observations and simulations. After each of the five eruptions, there was a cooling trend in the following three years (Fig. 1a). However, a global cooling pause against the monotonous cooling trend in global mean SAT also existed in GISTEMP during the first boreal winter after each eruption (hereafter the first boreal winter). Similar results can also be found in the Hadley Centre Climatic Research Unit surface temperature version 4 (HadCRUT4; Morice et al. 2012) and 20CRv2c. In this first boreal winter, the averaged global mean SAT recovery, which was identified by the difference between the original and 3-yr low-pass SAT, can reach 0.1 K for the average of the five eruptions. Figure 1b shows the MME mean of simulated global mean SAT anomalies for these tropical eruptions in the historical and AMIP runs of CMIP5 models. A significant cooling trend in the three years following each of these large eruptions was simulated in both the historical and AMIP runs (Fig. 1b), which is consistent with observations. Maximum cooling in the MME mean of the 32 historical models reaches 0.3 K in the boreal summer of the second posteruption year. The first-boreal-winter global cooling pause is reproduced by the MME mean of the AMIP runs, but is not adequately simulated by the MME mean of the historical models, which only exhibits a monotonous cooling trend. In summary, the MME mean of the CMIP5 models replicates global cooling in the three years following the large tropical eruptions but does not adequately simulate the first-boreal-winter response.
Figure 2 shows composite global SAT anomalies during the first boreal winter after the eruption in reanalysis, observations, and simulations. In 20CRv2c (Fig. 2a), a pause in global cooling during the first boreal winter appears as significant El Niño warming in the tropics and strong warming in Eurasia, which is consistent with previous finding (Robock and Mao 1992). Eurasian warming occurs during the first posteruption boreal winter after four of the five large tropical eruptions, with the 1963 Agung eruption being the only exception. Similar patterns can also be found in the observations of GISTEMP (Fig. 2b) and HadCRUT4 (not shown).
In the MME mean of the AMIP runs (Fig. 2c), the significant SAT response over the ocean is similar to that of the observations due to specified SST; a positive Pacific–North American (PNA)-like pattern in SAT is also seen from the tropical Pacific to North America, with positive–negative–positive–negative anomalies over the central Pacific, North Pacific, Canada, and southeastern United States. The Eurasian warming response in the first boreal winter, however, is much weaker in the models than in observations. Without the specified SST, the MME mean of the historical runs exhibited significant global cooling over low-latitude regions, with weak warming only occurring over the NH polar region (Fig. 2d). The El Niño–like response is not simulated by the MME mean of the historical runs.
In summary, global cooling in the three years following each of the five large tropical eruptions is simulated by both the historical and AMIP runs. The MME mean of the AMIP runs reproduces the first-boreal-winter SAT responses well, while that of the historical runs does not. One major difference between the historical and AMIP runs is the El Niño signal in the tropical Pacific, which is closely connected to global temperature anomalies through teleconnection caused by latent heating (Trenberth 2002). The Niño-3.4 index for these five tropical eruptions is shown in Fig. 3 for both observations and historical runs. In the three observational datasets, most of the eruptions were followed by El Niño events in the first boreal winter (Fig. 3a). The 1883 Krakatau eruption, being the only exception, was followed by a neutral ENSO response. Four of the five eruptions showed both global cooling pause and positive Niño-3.4 anomalies in the first boreal winter, while the 1883 Krakatau eruption only yielded a global cooling pause and a neutral ENSO response. The global cooling pause and El Niño response after the 1982 El Chichón eruption were stronger than those of the others (Figs. 1a and 3a). In the AMIP runs from 1979 to 2005, observational El Niño SST was specified as part of the forcing for the 1982 and 1991 boreal winters, while the MME mean of the historical runs from year (0) to year (2) only simulates a neutral ENSO response (Fig. 3b).
Since there is a close relationship between global mean SAT and ENSO (Trenberth 2002; Trenberth and Fasullo 2013), the observed and simulated ENSO and global mean SAT responses in the first boreal winter after each of these eruptions are exhibited in Fig. 4. Global mean SAT anomalies in the first boreal winter, accompanied by warm ENSO events, are nearly neutral or positive against the cooling trend in both observational mean and MME mean of the AMIP runs. Besides, the simulated global mean SAT anomalies in the MME mean of the AMIP runs are cooler than the observational mean for the 1982 El Chichón and 1991 Pinatubo eruptions although the SST field was specified. In the historical runs, this close relationship between global mean SAT and ENSO is also simulated, and the global mean SAT increases with the Niño-3.4 index at a ratio of 0.31, which is significant at the 95% confidence level. A few ensemble members, however, simulate the global cooling pause. The MME mean of the historical runs only replicates the neutral mode of ENSO and the negative global mean SAT anomalies after the eruptions, which is significantly cooler than the observed El Niño and neutral global mean SAT anomalies. Even though some models with weak El Niño–like responses, such as GFDL-ESM2M and ACCESS1.3, simulate weaker global cooling than the others do, none captures the observed global cooling pause during the first boreal winter (Fig. 4).
The MME mean of the historical runs fails to show the El Niño–like response in the first boreal winter after the tropical eruptions (Fig. 2d). It has been found that the initial ocean condition (IOC), namely, the ENSO state in the pre-eruption winter, is important to the volcano–ENSO relationship (Ohba et al. 2013; Pausata et al. 2016; Predybaylo et al. 2017). No significant El Niño–like response can be found in the first boreal winter after an eruption if there is already an El Niño in the boreal winter before the eruption (Liu et al. 2018a). The Niño-3.4 index in the last boreal winter before the eruption was used as the IOC. To isolate El Niño from volcanic forcing, relative SST (Khodri et al. 2017) is employed here. In observations, the IOC of each eruption happened to be a La Niña or neutral state (Fig. 3a), the composite El Niño during the first boreal winter after the eruption is developed from a weak La Niña IOC, and an El Niño signal with westerly wind anomalies appears in the boreal summer and peaked in winter (Fig. 5a).
In the MME mean of the historical runs, there is no significant warming signal in the eastern equatorial Pacific in the first boreal winter after the eruptions (Fig. 5b). When the IOC is not in an El Niño state (Fig. 5c), the MME mean of the historical runs simulates the first-boreal-winter El Niño–like response that is significant in the central Pacific. However, the simulated response is much weaker than that in the observations, suggesting relatively weak air–sea interaction in these models. For an El Niño IOC (Fig. 5d), the response is determined by the evolution of El Niño, and a La Niña–like response is simulated in the first boreal winter, which is consistent with previous finding (Liu et al. 2018a; Predybaylo et al. 2017). Similar results can also be seen from large ensemble simulations of the CESM (Figs. 5f–h). The El Niño–like response in the first boreal winter after the eruptions is simulated with a La Niña or neutral IOC. In the CESM-LE, the El Niño–like response peaks in the second year after the eruptions.
In summary, the MME mean of the historical runs reproduces the El Niño–like response in the first boreal winter after the tropical eruptions when the IOC is not in an El Niño state, which is consistent with observations. The simulated El Niño–like response, however, is much weaker than that in the observations.
4. Eurasian SAT and polar vortex responses
In the first boreal winter after each of the five tropical eruptions, a Eurasian warming occurred, which contributed to the pause in the 3-yr global cooling trend following the eruption (Fig. 2b). This Eurasian warming was also noted in previous works based on observations (Robock 2002; Robock and Mao 1992) and reconstructions (Fischer et al. 2007). The Eurasian SAT has been found to be affected by the polar vortex, namely, a strong stratospheric polar vortex related to the NAO can induce significant Eurasian warming through planetary wave reflection and jet stream disturbance (Butler et al. 2014; Perlwitz and Graf 1995; Thompson and Wallace 2001). As shown in Fig. 6, the high correlation between Eurasian SAT and stratospheric polar vortex anomalies is found in both reanalysis datasets and the historical and AMIP runs of the CMIP5 models, which indicates that the CMIP5 models can simulate the polar vortex–Eurasian SAT relationship well.
Although the MME mean of the AMIP runs shows the first-boreal-winter global cooling pause after the eruption, this simulated global mean SAT is colder than that in observations (Fig. 4). The main reason is that the simulated Eurasian warming response in the AMIP runs is weaker than the observed one (Figs. 2b,c). The historical runs also fail to simulate the strong Eurasian warming response (Fig. 2d). The failure in reproducing the Eurasian warming response may result from poor simulation of stratospheric response, especially of the polar vortex (Driscoll et al. 2012).
Since the reanalysis of 20CRv2c only has assimilated surface observations (Compo et al. 2011), the reliability of its stratospheric variables is questionable. The reanalysis of polar vortex index is obtained from JRA-55 that starts in 1958. In the first boreal winter following each of the three eruptions after 1958, the observed Eurasian SAT anomaly was positive and was associated with an enhanced polar vortex (Fig. 7), which is consistent with previous findings (Zambri and Robock 2016). In both the historical and AMIP runs, the Eurasian SAT anomaly is significantly correlated to the anomalous polar vortex, with the correlations being 0.47 and 0.50, respectively; both are significant at the 95% confidence level (Figs. 7a,b). In the MME mean of the historical runs, the simulated Eurasian warming and polar vortex enhancement in the first boreal winter are significantly weaker compared to the observed ones (Fig. 7a). The strong polar vortex, exceeding 0.5 standard deviations of the reanalysis in the period of 1958–2005, is only simulated in HadGEM2-AO and the three NCAR models (CESM1-CAM5-fV2, CESM1-WACCM, and CESM1-BGC). Although each of these models has only one ensemble member, it may be good at representing stratospheric processes; the reason for this requires further study.
In the AMIP runs (Fig. 7b), which are forced by specified SST, only NorESM1-M, GFDL-HIRAM-C360, and two BCC models (BCC-CSM1.1-m and BCC-CSM1.1) simulate the strong polar vortex enhancement and Eurasian warming. It can be seen from the MME mean of the AMIP runs that the Eurasian warming and polar vortex enhancement are relatively weaker than those in the observations, which means that most models that participated in the AMIP runs underestimate the direct volcanic aerosol effect in enhancing the polar vortex.
Figure 8 shows the composite distribution of polar vortex responses in different reanalysis datasets and in CMIP5 simulations. In JRA-55 (Fig. 8a), a strong positive polar vortex anomaly was observed at 50 hPa during the first boreal winter after each of the three eruptions since 1958. The center of the polar vortex anomaly was over the Arctic region. Due to the tropical warming and polar cooling in the stratosphere, strong westerly wind anomalies appeared near 70°N (Fig. 8b). All the three eruptions, namely, the Agung, El Chichón, and Pinatubo eruptions, were associated with similar polar vortices, but the polar vortex of the Pinatubo eruption was weakly increased. This strengthened polar vortex can also be seen in the ERA-40 reanalysis after the 1982 El Chichón and 1991 Pinatubo eruptions (Zambri and Robock 2016) and in the NCEP reanalysis during the two posteruption winters (Graf et al. 2014). Compared to the result based on JRA-55, the positive polar vortex anomaly in the first boreal winter after the eruption in 20CRv2c is very weak, and no significant anomalous geopotential height and temperature signals can be detected over the polar stratosphere (Figs. 8c,d). Without assimilation of upper-level variables, the NCEP Global Forecast System model used for 20CRv2c still has difficulty in capturing observed volcano-induced polar vortex change.
In the MME mean of the AMIP runs (Fig. 8e), the polar vortex anomaly is very weak, and a dipole pattern after the eruption is simulated with significant positive geopotential height anomaly over Alaska and weak negative anomaly over the Greenland Sea, which is not consistent with the results based on JRA-55 (Fig. 8a). The historical runs suffer the same problem (Fig. 8g). Although the warming signal in the tropical stratosphere is simulated, the polar vortex is not enhanced significantly and the westerly anomalies around the North Pole are much weaker than those based on JRA-55 (Figs. 8f,h), which means the direct volcanic aerosol effect on enhancing the polar vortex is underestimated in the CMIP5 models. The strength of volcanic forcing should be large enough to detect the significant polar vortex response (Bittner et al. 2016). Our results based on the strong 1883 Krakatau and 1991 Pinatubo eruptions also showed that the response of the polar vortex is underestimated in the historical runs of the CMIP5 models (not shown).
For all of these 32 CMIP5 models, when a strong response is defined by exceeding 0.5 standard deviations of observations or reanalysis, none of them, either with a La Niña or neutral IOC, simulate the observed El Niño and strong positive polar vortex anomalies at the same time in the first boreal winter after the tropical eruption (Fig. 9). Among the 97 ensemble members of these 32 models, only four ensemble members simulate El Niño and strong positive polar vortex anomalies simultaneously.
In summary, the stratospheric polar vortex is mainly affected by direct volcanic aerosol-induced stratospheric warming in the reanalysis. Despite being reproduced by both the historical and AMIP runs of the CMIP5 models, the enhanced polar vortex in the first boreal winter after the eruption is very weak compared to that based on the reanalysis. Not only the CMIP5 models but also the NCEP Global Forecast System model underestimates the direct volcanic aerosol effect, and neither of them adequately captures the role of volcanic forcing in changing the polar vortex, which are consistent with previous findings (Marshall et al. 2009). Because a higher stratospheric resolution does not increase a model’s ability to simulate the response to volcanic forcing (Charlton-Perez et al. 2013; Marshall et al. 2009), deficient nonlinear stratosphere–troposphere interaction (Stenchikov et al. 2006) may explain these models’ limitation.
5. Discussion on internal variability–induced global cooling pause
We cannot conclude that this global cooling pause is totally excited by tropical eruptions only based on these five samples, because internal variability may be a potential cause (Polvani et al. 2019). Forty-one ensemble members of the CESM-LE are thus analyzed for the latest three eruptions since 1920 to address this question.
The significant linear relationship between global mean SAT and Niño-3.4 index anomalies and that between the polar vortex and Eurasian SAT anomalies are simulated by the 41 ensemble members of the CESM-LE, with correlation coefficients being 0.76 and 0.54, respectively; both are significant at the 95% confidence level (Fig. 10). The simulated global mean SAT and Niño-3.4 index anomalies in the ensemble mean of the CESM-LE in the first boreal winter are −0.13 and 0.12 K, respectively (Fig. 10a), which are much colder than the observations, although they are warmer than those simulated by the CMIP5 historical runs (Fig. 4). The simulated polar vortex enhancement and Eurasian warming in the ensemble mean are also much weaker than the observed ones (Fig. 10b).
Among these 41 ensemble members of the CESM-LE, four simulate significant global cooling pause as well as El Niño responses (Fig. 10a), and other four replicate significant polar vortex enhancement and Eurasian warming (Fig. 10b), which exceed 0.5 standard deviation of observations or reanalysis. For the ensemble mean of the CESM-LE (Fig. 11a), tropical cooling is reproduced, but the central equatorial Pacific warming is weak. Significant warming is only simulated over the Arctic region (Fig. 11a), and the polar vortex enhancement is weak (Fig. 11b). Ensemble members that simulate the global cooling pause reproduce SAT responses similar to those in observations except that the Eurasian warming is shifted northward (Fig. 11c), but they fail to simulate strong polar vortex enhancement (Fig. 11d). For the ensemble members that capture strong polar vortex enhancement and Eurasian warming (Figs. 11e and 8f), the global cooling pause, however, is not simulated, especially over the tropics (Fig. 11e).
The ensemble mean of the CESM-LE fails to show the global cooling pause; also, the simulated polar vortex enhancement and Eurasian warming are very weak. Although some ensemble members capture the global cooling pause or polar vortex enhancement, they do not simulate both well.
Whether the global cooling pause in the first boreal winter is a response to a tropical eruption or is internal variability is still an open question. Strong polar vortex enhancement in reanalysis data is argued to result from internal variability (Bittner et al. 2016; Polvani et al. 2019). We conclude that the global cooling pause is a response to the tropical eruption based on long-term reconstruction evidence that shows a significant El Niño response (Adams et al. 2003; Liu et al. 2018a,b) and a significant overall warm anomaly over northern Europe (Fischer et al. 2007) in the first boreal winter after the tropical eruption, and based on physical processes reviewed in the introduction, including the dynamic thermostat mechanism and the temperature gradient mechanism.
In this study, we investigated global SAT responses during the first boreal winter after the five major tropical volcanic eruptions since 1870 compared simulations of CMIP5 historical and AMIP runs and of CESM-LE to observations and reanalysis. In observations, the global cooling pause, which is due to tropical El Niño–like warming and Eurasian warming, is seen during the first posteruption boreal winter, thereby forming a deviation from the 3-yr global cooling trend after the tropical eruptions. Direct volcanic aerosol-induced warming in the tropical lower stratosphere and cooling in the polar stratosphere increase the meridional temperature gradient and enhance the stratospheric polar vortex, resulting in strong Eurasian surface warming (Figs. 8a,b).
In the AMIP runs, the first-boreal-winter global cooling pause is simulated with specified SST forcing. In the historical runs, a monotonous global cooling trend rather than a global cooling pause is produced. The IOC is one key to simulate the first-boreal-winter El Niño response, as the latter is reproduced when the pre-eruption IOC is not in an El Niño state, but tropical cooling is still produced since the simulated El Niño response is very weak compared to that in the observations. Although the CMIP5 models show the Eurasian SAT response to the polar vortex (Fig. 6), the simulated Eurasian warming and positive stratospheric polar vortex responses to the eruptions are very weak in both AMIP and historical runs (Fig. 8), which means that the CMIP5 models underestimate the direct volcanic aerosol effect on the stratosphere. The majority of the CMIP5 coupled models cannot simulate both El Niño and enhanced polar vortex responses (Fig. 9). The abovementioned processes controlling the global cooling pause in the first posteruption boreal winter are summarized in Fig. 12. Although this global cooling pause is simulated by some ensemble members of the CESM-LE as an internal mode (Figs. 10 and 11), we conclude that it is a response to tropical volcanic eruption according to the reconstruction analysis.
Different from the observations, the models fail to simulate the first-boreal-winter El Niño after a tropical eruption (Fig. 5). This suggests that the CMIP5 models underestimate the role of the thermostat mechanism (Clement et al. 1996) or the Bjerknes feedback (Bjerknes 1969) in exciting El Niño response. Neither the CMIP5 models nor the NCEP Global Forecast System model used for 20CRv2c can simulate the strong polar vortex in the first boreal winter after a tropical eruption. All these deficiencies in simulating global SAT responses to the tropical eruptions call for urgent improvement of models; and the upcoming model intercomparison project on the climate response to volcanic forcing (VolMIP; Zanchettin et al. 2016) offers a good opportunity to discuss these deficiencies. Analysis of climate responses to large volcanic eruptions provides an additional perspective to evaluate our models’ sensitivity to external forcing.
This work is supported by the Natural Science Foundation of China (41975107), the Public Science and Technology Research Funds Projects of Ocean (201505013), the National Science Foundation of the United States (AGS-1540783), the National Natural Science Foundation of China (41971108), and the Strategic Priority Research Program of Chinese Academy of Sciences (XDA20060401). We acknowledge the World Climate Research Programme’s Working Group on Coupled Modelling, which is responsible for CMIP, and we thank the climate modeling groups for producing and making available their model outputs. This paper is ESMC Contribution No. 297.
Denotes content that is immediately available upon publication as open access.