3.1 Global mean volcanic response

Figure 2 shows the global mean net top-of-atmosphere (TOA) radiation flux response (DCPP-A minus DCPP-C) to the individual volcanic eruptions and the mean of the three. The TOA flux is computed as the incoming shortwave (i.e. the rsdt CMIP6 variable) minus the outgoing shortwave (i.e. rsut) and longwave radiation (i.e. rlut) fluxes. The climate predictions show a statistically significant post-volcanic decrease in global mean TOA which is consistent across the models, generally peaking in the first boreal winter. There are differences, however, in the magnitude of the response and the recovery following the three volcanic eruptions, coherent with the volcanic aerosol forcing (Fig. 1). After the eruption of Agung, the negative TOA anomalies reach approximately −1.7 ± 0.4 W m⁻² (model mean ± the inter-model standard deviation) and recover within the next months. The TOA response for the eruption of El Chichón is the weakest and the only one in which the model spread overlaps with the zero line, with the ensemble mean response reaching approximately −0.9 ± 0.2 W m⁻² and recovering within the next year. The response to the eruption of Pinatubo shows the strongest TOA negative anomalies of the three, reaching approximately −2.3 ± 0.3 W m⁻² and recovering in approximately 2 years. These results are consistent with those of Zanchettin et al. (2022).

Figure 2Global mean top-of-atmosphere radiation response (W m⁻²) to the volcanic eruptions (DCPP-A minus DCPP-C). The ensemble mean for each model and the multi-model mean are shown. The shading is the multi-model member spread calculated as the 10th and 90th percentiles of the entire ensemble. Filled squares at the bottom part of the figure indicate statistically significant differences (see the Methods section). The vertical dashed grey lines indicate the approximate time of the eruptions. The data for CESM1-1-CAM5-CMIP5 were not included since they are not available. Figure inspired by Zanchettin et al. (2022).

Download

In response to the negative TOA anomalies following the volcanic eruptions, the global mean surface temperature decreases (Fig. 3). The cooling anomalies peak in the second year and subsequently recover in approximately 5 years. Similarly to the TOA response, the prediction systems generally show a comparable global mean surface temperature response for the individual eruptions, although noisier, and the differences are larger across the eruptions. For the eruptions of Agung and El Chichón, the minimum global mean surface temperature anomalies are comparable (the model mean ± the inter-model standard deviation is −0.21 ± 0.05 and −0.17 ± 0.04 °C, respectively) even though the TOA anomalies differences are larger, indicating there are potential non-linearities in the response or multiple mechanisms at play. The eruption of Mount Pinatubo shows the largest global mean surface temperature anomaly of −0.35 ± 0.06 °C (model mean ± the inter-model standard deviation), coherent with the larger TOA anomalies. There is considerable inter-model spread in the response to both the Agung and El Chichón eruptions, while for Pinatubo the models show stronger agreement, probably because of the greater signal-to-noise ratio. The good overlap between the temperature anomalies reported here and found in Zanchettin et al. (2022) (multi-model mean of about −0.33 °C) for the eruption of Pinatubo suggests that the influence from the initial conditions and the background state is small, at least for the global mean.

Figure 3Global mean surface air temperature response (°C) to the volcanic eruptions (DCPP-A minus DCPP-C). The ensemble mean for each model and the multi-model mean are shown. The shading is the multi-model member spread calculated as the 10th and 90th percentiles of the entire ensemble. Filled squares on the bottom part of the figure indicate statistically significant differences (see the Methods section). The vertical dashed grey lines indicate the approximate time of the eruptions.

Download

Figure 4 shows the changes in ocean heat content (OHC) in the upper 300 m, an integrated variable with large inertia for which the volcanic signals can persist longer in time (e.g. Stenchikov et al., 2009). The volcanic signals are detected throughout the whole decade following the eruptions, with minimum OHC values peaking around the fourth year following the eruption and a slow recovery that still yields statistically significant anomalies in the ensemble mean by the end of the forecasts in each of the individual eruptions. Likewise, all systems simulate long-lasting significant responses, with some inter-model differences in the peak timing and the recovery rate.

Figure 4Global ocean heat content response (J m⁻²) to the volcanic eruptions (DCPP-A minus DCPP-C). The ensemble mean for each model and the multi-model mean are shown. The shading is the multi-model member spread calculated as the 10th and 90th percentiles of entire ensemble. Filled squares at the bottom part of the figure indicate statistically significant differences (see the Methods section). The vertical dashed grey lines indicate the approximate time of the eruptions.

Download

Having characterised the global mean response to the three volcanic eruptions, we compare the global mean near-surface air temperature anomalies predicted in both the DCPP-A and DCPP-C hindcasts against observations. Figure 5 shows that overall the DCPP-A predictions, which include the volcanic forcing, reproduce the HadCRUT5 temperature anomalies more closely than the DCPP-C. The DCPP-A hindcasts generally simulate the observed cooling tendencies in the initial years following the eruptions and the observed anomalies (from HadCRUT5) are generally within the 10th and 90th percentiles. In contrast, the DCPP-C hindcasts are generally warmer and observations tend to fall more frequently below the 10th and 90th percentiles. This is shown by lower RMSE of the DCPP-A hindcasts in comparison to the DCPP-C (Table S1). In particular, the volcanic forcing following the eruption of Pinatubo is remarkably important to simulate the observed global mean surface temperature variability in the early 90s, as without the forcing the observations are out of the models' envelope in the 2 years following the eruption.

Figure 5Monthly mean global near-surface temperature anomalies (°C) of the predictions initialised in 1962, 1981 and 1990 for the DCPP-A (with volcanic forcing) and DCPP-C (without volcanic forcing) experiments. HadCRUT5 is used as the observational reference (dashed line). The anomalies have been computed with respect to the period 1970–2005 (see the Methods section for further information). The shading is the multi-model member spread calculated as the 10th and 90th percentiles of the entire ensemble.

Download

3.2 Spatiotemporal characteristics of the volcanic response

The surface air temperature response in the first year (computed June–May to characterise the post-eruption anomalies) following the volcanic eruptions (Fig. 6) shows a distinct pattern, largely consistent across the volcanic eruptions (as shown by the shading in Fig. 6). It is characterised by a generalised cooling that is largest in the tropics and subtropics and a warming in the Eurasian Arctic sector (discussed further in Sect. 3.3). The magnitudes of the anomalies tend to be greater over land than over the ocean, in agreement with response to radiative forcing (e.g. Eyring et al., 2021). Despite the overall similarities in the surface temperature response among the eruptions, there are some differences which can be explained by the magnitude and geographical pattern of the TOA anomalies following each eruption (Fig. S1). For the eruption of Agung the TOA anomalies mainly occur in the Southern Hemisphere (Fig. S1a), and therefore the surface temperature anomalies are larger over the tropics and Southern Hemisphere (Fig. 6a). In contrast, for the eruption of El Chichón the TOA anomalies occur mainly in the Northern Hemisphere (Fig. S1b) and the surface temperature anomalies are larger over the tropics and the Northern Hemisphere (Fig. 6b). The eruption of Pinatubo caused a hemispherically symmetric response (Fig. S1c) with surface temperature anomalies that were stronger over the tropics and at high latitudes (Fig. 6c).

Figure 6Model mean near-surface air temperature (°C) response (DCPP-A minus DCPP-C) during the first year following the volcanic eruptions (June–May). Hatching indicates statistically significant anomalies, while the shading indicates model agreement (see the Methods section).

For years 2–5 (computed as June–May), the cooling spreads worldwide in response to the volcanic forcing (Fig. 7a–d). There are evident differences among the eruptions, associated with the magnitude and geographical pattern of the forcing as described previously. For the eruption of Agung the largest surface temperature anomalies are located over the tropics and especially in the Southern Hemisphere sector (Fig. 7a), while following the eruption of El Chichón the strongest surface temperature anomalies are located over the polar latitudes of the Northern Hemisphere, probably linked to the Arctic amplification phenomenon (Fig. 7b). The eruption of Pinatubo shows stronger surface temperature anomalies worldwide, albeit with maxima over the Arctic like for El Chichón (Fig. 7c). The mean response to the three eruptions is consistent with the one to Pinatubo, with surface temperature anomalies that are statistically significant worldwide (Fig. 7d).

At later forecast times (years 6–9) the surface temperature anomalies largely decrease in magnitude (Fig. 7e–h), consistent with the recovery timescale after the eruptions (Fig. 3). The response on these timescales is partly related to the magnitude of the eruption. For the eruption of Agung we find that negative surface temperature anomalies persist over the tropics and the Southern Hemisphere (Fig. 7e). For the eruption of El Chichón, which is the weakest of the three, there are barely any significant anomalies remaining (Fig. 7f). Finally, for the eruption of Pinatubo, we find some regions with significant anomalies which are particularly strong over the Arctic (Fig. 7g).

Figure 7Multi-model mean near-surface air temperature (°C) response (DCPP-A minus DCPP-C) during years 2–5 (first row) and years 6–9 (second row) following the eruptions of Mount Agung (1963), El Chichón (1982) and Mount Pinatubo (1991) as well as the mean of the three volcanoes (left to right). The annual mean is defined from June to May. Hatching indicates statistically significant anomalies, while shading indicates model sign consistency (see the Methods section).

To determine the importance of volcanic forcing at the regional scale, we evaluate the predicted surface temperature anomalies in the three DCPP-A and DCPP-C hindcast sets against the HadCRUT5 observational dataset. We focus on the multi-model ensemble mean for year 1, years 2–5 and years 6–9. Overall, we find that the volcanic forcing only leads to a generalised improvement in forecast years 2–5, as indicated by the lower global RMSE the three DCPP-A hindcasts with respect to DCPP-C (Table S2; by 12 % for Agung and El Chichón and 16 % for Pinatubo). Figure 8 shows that the multi-model mean surface temperature anomaly patterns for years 2–5 (when the volcanic radiative impact is strongest) predicted by the DCPP-A hindcasts are cooler than the DCPP-C, as expected, and closer to the observations. In comparison with the observed anomaly patterns, the multi-model predicted pattern is smoother and does not reproduce most of the regional variations. This is probably because at this forecast range the multi-model pattern is mostly due to the radiatively forced response and the observed regional variations are due to noise. For the first forecast year and forecast years 6–9, although the volcanic forcing has a distinct surface temperature imprint following the eruptions (Figs. 6 and 7), the fact that we do not find a detectable improvement in the DCPP-A hindcasts over DCPP-C (Table S2) might be because the local volcanic response is overwhelmed by internally generated variability and/or the regional response to the volcanic forcing is not correctly simulated by the models.

Figure 8Near-surface air temperature anomalies (°C) for years 2 to 5 of the multi-model mean predictions initialised in 1962, 1981 and 1990 for the DCPP-A (with volcanic forcing) and DCPP-C (without volcanic forcing) experiments and HadCRUT5. The anomalies have been computed with respect to the period 1970–2005 (see the Methods section). Note that the anomalies are computed from June–May.

3.3 Northern Hemisphere atmospheric response

The models consistently simulate a surface warming over northern Eurasia and/or the Barents–Kara seas in the first year after the eruptions (Fig. 6), although it is only statistically significant for Agung and the composite of the three eruptions. A focus on the seasonal surface temperature changes reveals that after the eruptions of Agung and Pinatubo, the surface warming over northern Eurasia and the Barents–Kara seas reaches its maximum magnitude and extent in the first winter (Fig. 9a and c), although only small regions show statistically significant anomalies. Interestingly, for the Agung eruption the warming is already present in the first JJA and SON (Fig. S2a–d). In contrast, following the eruption of El Chichón there are no warm anomalies present over the Eurasian continent in the first DJF (Fig. 9b), which shows cool anomalies instead (not statistically significant) and very weak and localised warm anomalies over the Barents–Kara seas. Due to the disparate responses among the eruptions, the composite of the three eruptions shows positive yet not significant anomalies in both DJF and the earlier seasons (Fig. S2). These results are consistent with studies suggesting that the Eurasian warming might be too weak and hence overwhelmed by unforced variability to all major eruptions (DallaSanta and Polvani, 2022).

Coherent with the surface temperature anomalies, for the Agung and the Pinatubo eruptions the sea level pressure anomalies feature a positive NAM-like pattern during the first post-eruption winter (Fig. 9e and g), not occurring for the eruption of El Chichón (Fig. 9f). The positive NAM-like pattern in the first post-eruption winter is associated with the warming of the tropical lower stratosphere in the months following the volcanic eruptions (Fig. 10a–d), which increases the poleward temperature gradient and might explain the acceleration of the polar vortex (Fig. 10e–h), in line with previous studies (e.g. Hermanson et al., 2020). The lower-stratospheric temperatures in the tropics show a strong increase following the three eruptions, peaking approximately 6 months after the eruption and lasting for 3 years (Fig. 10a–d), thereby increasing the poleward temperature gradient (not shown). Four of the forecast systems (CanESM5, CMCC-CM2-SR5, EC-Earth3 and HadGEM3-GC31-MM) cluster together, while IPSL-CM6A-LR and CESM1-1-CAM5-CMIP5 tend to simulate greater anomalies in response to the eruptions. In the case of CESM1-1-CAM5-CMIP5 this could be because this model was forced with CESM-specific volcanic forcing (Ammann et al., 2003) rather than CMIP6. The temperature anomalies are coherent with the magnitude of the volcanic forcings; the multi-model mean reaches 4.7 ± 0.7, 2.4 ± 1.4 and 5.6 ± 1.7 °C for the eruptions of Agung, El Chichón and Pinatubo, respectively.

Accompanying the stratospheric warming there is a detectable acceleration of the polar vortex following the eruptions of Agung and Pinatubo in the first winter in both the model mean and the eruption mean (Fig. 10e, g and h). Not all individual models show this acceleration though, probably because a larger ensemble size is needed to overcome the low signal-to-noise ratio and detect the response. Furthermore, a stronger stratospheric temperature response does not necessarily produce a stronger polar vortex response. In the case of the eruption of El Chichón the warm anomalies in the tropical lower stratosphere are weaker (Fig. 10b), which does not seem to result in an acceleration of the polar vortex and explains why a positive phase of the NAM and its associated warming over Eurasia are not simulated (Fig. 10f). The sea level pressure anomalies do resemble a positive NAM-like pattern in autumn, when weak warm anomalies are present over northern Eurasia (Figs. S2f and S3f), but this is not accompanied by a detectable acceleration of the polar vortex.

Figure 9Multi-model and multi-eruption response (DCPP-A minus DCPP-C) of surface air temperature (a–d) and sea level pressure (e–h) in the first boreal winter (DJF) following the volcanic eruptions. The ensemble mean for each model and the multi-model mean are shown. Shading indicates statistically significant anomalies.

Figure 10Stratospheric air temperature in the tropics (30° N–30° S at 50 hPa) and polar vortex (average zonal velocity over 55–75° N at 50 hPa) response (DCPP-A minus DCPP-C) following the volcanic eruptions. The ensemble mean for each model and the multi-model mean are shown. The shading is the multi-model member spread calculated as the 10th and 90th percentiles. Filled squares on the bottom part of the figure indicate statistically significant differences (see the Methods section).

Download

3.4 Response in the Pacific Ocean

Many studies have reported the impact of volcanic eruptions on the variability of the Pacific Ocean from seasonal to decadal timescales. We start by documenting the impact of volcanic eruptions on ENSO in our simulations. To isolate the dynamical response of ENSO from the surface cooling effect we define the Niño 3.4 SST index relative to the tropical SST mean (20° N–20° S) as in Khodri et al. (2017). Figure 11 shows that there is a large uncertainty in the relative Niño3.4 response to the volcanic eruptions and no consistent response across the individual models. This is partly due to the small ensemble size for the individual models (10 ensemble members), since large ensembles have been previously shown to be required to detect such signals (e.g. Ménégoz et al., 2018a; Hermanson et al., 2020). Only when considering the model mean and multi-eruption composite do we find a clear and statistically significant response. It is characterised by the development of weak El Niño-like conditions in the year of the eruptions which peak in the following year and then a transition to La Niña-like conditions in the second and third years following the eruption. There are, however, some differences among the volcanic eruptions, as the multi-model El Niño response is stronger and significant only for the eruptions of Agung and Pinatubo (Fig. 11a and c), while the delayed La Niña response is clear and significant for Pinatubo and marginally significant (i.e. only during a season) for El Chichón (Fig. 11b). From the individual models CMCC-CM2-SR5 is the only one consistently showing the El Niño-like and La Niña-like responses after the three eruptions; all the rest are not significant in almost all cases. HadGEM3-GC31-MM and CESM1-1-CAM5-CMIP5 also show a significant positive ENSO response in the multi-eruption composite, with CESM1-1-CAM5-CMIP5 showing a subsequent significant La Niña-like response. None of the other three systems simulate the El Niño/La Niña response. Note that the results for the IPSL-CM6A-LR model contrast with those of Zanchettin et al. (2022), who show in their idealised Pinatubo experiments (volc-pinatubo-full experiment) that this model does simulate the El Niño/La Niña responses.

Figure 11Relative Niño 3.4 index response following the eruptions of (a) Mount Agung (1963), (b) El Chichón (1982) and (c) Mount Pinatubo (1991) as well as (d) the mean of the three eruptions. Filled squares on the bottom part of the figure indicate statistically significant differences (see the Methods section). The ensemble mean for each model and the multi-model mean are shown. The shading is the multi-model member spread calculated as the 10th and 90th percentiles. The vertical dashed lines indicate the approximate time of the eruptions.

Download

To further explore the development of the El Niño-like conditions we look at the model mean composite responses in surface air temperature and precipitation in the tropics in JJA and DJF of the first and second years following the eruptions. The eruption mean response shows the development of a Niño-like state from the first JJA to the second DJF (Fig. 12), comparable to the results shown in Hermanson et al. (2020). In the year of the eruption, the tropics cool, especially over the continents, and there is a reduction in precipitation over Africa and the Maritime Continent (Fig. 12b and d). In the following year's winter, the eastern and central equatorial Pacific show enhanced warming relative to the rest of the tropics, accompanied by the characteristic El Niño-like precipitation pattern (Fig. 12g and h).

While combining the different eruptions is a way to increase the ensemble size and improve the detection of weak signals, in this case it is probable that due to the different characteristics of each eruption, timing and background climate state, the mechanisms at play are not the same for each eruption and therefore their impact on ENSO also changes. It is beyond the scope of this study to determine which mechanisms might dominate the response in the individual eruptions since these seem to vary from one model to another (Fig. 11a-c). Nonetheless we find that overall, for the eruptions of Agung and Pinatubo, the El Niño-like state develops and peaks in the first year following the eruptions, while for the eruption of El Chichón the El Niño-like state occurs in the same year of the eruption. In the case of the Agung, initially a weak La Niña-like conditions develop in the year of the eruption (Fig. S4a and b) accompanied by a northern shift of the Intertropical Convergence Zone (ITCZ) (Fig. S5a and b); then the conditions shift to El Niño-like with a lack of cooling relative to the rest of the tropics (Fig. S4c and d) and accompanied by increased precipitation (Fig. S5c and d). For the eruption of El Chichón, we find an El Niño-like state developing and peaking in the year of the eruption (Fig. S4e and f). The responses to the eruptions of Agung and El Chichón are broadly consistent with the results of Pausata et al. (2020), who show using idealised simulations with NorESM1-M that the ENSO response to NH and SH eruptions differs due to shifts in the ITCZ as a result of the asymmetric hemispheric cooling. Finally, for the eruption of Pinatubo, which due to its intensity and location had an approximately hemispherically symmetric forcing and also induced the strongest cooling anomalies (Fig. S4j–m), El Niño-like conditions develop and peak in the year following the eruption. In this case the response is in broad agreement with the mechanism proposed in Khodri et al. (2017), which starts with reduced precipitation over tropical Africa (Fig. S5j and k), leading to the propagation of anomalous atmospheric Kelvin waves that weaken the trade winds over the western Pacific.

The ocean state has been found to be another relevant factor modulating the ENSO response in climate model simulations (e.g. Predybaylo et al., 2020). Observations show that in the months prior to the three volcanic eruptions considered in this study, El Niño phases were already developing, peaking later in the first winter following the eruptions (e.g. McGregor et al., 2020). Figure S6 shows that the observed ENSO anomalies in the first months after the eruptions were reasonably well predicted for Agung and less so for El Chichón, but not for Pinatubo, in both DCPP-A and DCPP-C sets. The very small differences identified between the two forecast ensembles suggest that the volcanic forcing has a weak impact, which could either imply that the observed ENSO signal in the first year was not forced by the eruptions or that the models systematically fail to simulate the true mechanism of response. At later forecast times, the models also misrepresent the observed ENSO anomalies.

Figure 12Multi-model and multi-eruption surface air temperature (°C) and precipitation (mm d⁻¹) responses following the eruptions in the tropics. Hatching indicates statistically significant anomalies (see the Methods section).

Previous studies have shown that including the volcanic forcing in decadal climate predictions can degrade the skill of SST over the central-eastern tropical Pacific Ocean on multiyear to decadal timescales (Timmreck et al., 2016; Wu et al., 2023). In these predictions, the tropical Pacific SST (averaged over 20° S–20° N, 120° E–80° W as in Wu et al., 2023) experiences a net cooling during the first 4 years following the eruptions, although with important differences among the models and volcanic eruptions regarding the magnitude of the cooling (Fig. S7). At later forecast times (years 6–10) there is no significant response to the volcanic forcing, with no evident differences between the ensembles. To evaluate whether the volcanic forcing has a detrimental impact on these predictions we compare the forecast anomalies of the DCPP-A and DCPP-C hindcasts with SST observations. Figure 13 shows that for the eruptions of Agung and El Chichón the differences between the DCPP-A and DCPP-C hindcasts are small (Table S3), although including the volcanic aerosols may improve the initial years of the forecasts. For the eruption of Agung, 1964 is too warm unless volcanic aerosols are included (Fig. 13a and b), and for the eruption of El Chichón, the multi-model tendency for 1983–1985 is better reproduced in DCPP-A compared to DCPP-C (Fig. 13c and d). In contrast, for the eruption of Pinatubo, including the volcanic forcing negatively impacts the prediction initialised in 1990 in all systems by causing a local cooling that was not observed. This is also evident from the RMSE values shown in table S3.

Figure 13Tropical Pacific (20° S–20° N, 160° E–80° W) SST anomalies (°C) in the predictions initialised in 1962, 1981 and 1990 for the DCPP-A (with volcanic forcing) and DCPP-C (without volcanic forcing) experiments. HadISSTv4 is used as the observational reference (dashed line). The anomalies have been computed with respect to the period 1970–2005 (see the Methods section for further information). The ensemble mean for each model and the multi-model mean are shown. The shading is the multi-model member spread calculated as the 10th and 90th percentiles of the entire ensemble.

Download

3.5 Response in the North Atlantic Ocean

The North Atlantic Ocean is a region where the impact of volcanic eruptions has been shown to persist on decadal timescales (e.g. Ortega et al., 2012; Swingedouw et al., 2015, 2017; Hermanson et al., 2020; Borchert et al., 2021; Fang et al., 2021). Volcanic eruptions will always tend to cause a direct cooling of North Atlantic SSTs via the worldwide reduction in incoming shortwave radiative fluxes (e.g. Swingedouw et al., 2017). As shown for the global mean surface temperature, in the North Atlantic Ocean, the SST decreases following the eruptions with a subsequent recovery (Fig. S8). There are important differences regarding the magnitude of the cooling and the recovery time, which are larger across the volcanic eruptions than across models. Comparing the predicted SST anomalies of the DCPP-A and DCPP-C hindcasts with SST observations shows that overall the volcanic forcing positively impacts the predictions (also shown by the RMSE values in Table S4). This is particularly relevant following the eruption of Pinatubo, for which including the volcanic forcing is necessary to realistically simulate the North Atlantic SST variability in all models, at least in the first few years (Fig. 14).

Figure 14North Atlantic (0–60° N, 80° W–0°) SST anomalies (°C) in the predictions initialised in 1962, 1981 and 1990 for the DCPP-A (with volcanic forcing) and DCPP-C (without volcanic forcing) experiments. HadISSTv4 is used as the observational reference (dashed line). The anomalies have been computed with respect to the period 1970–2005 (see the Methods section for further information). The ensemble mean for each model and the multi-model mean are shown. The shading is the multi-model member spread calculated as the 10th and 90th percentiles of all members.

Download

To isolate the changes that are specific to the North Atlantic Ocean, which might additionally arise from internal variability processes, we compute the AMV index as a standardised anomaly relative to the global SST mean (between 60° N–60° S) (Trenberth and Shea, 2006). We find that for this AMV index there is no significant impact due to the volcanic forcing (Fig. S9). As previous studies have shown, the AMOC plays a key role in modulating the AMV, associated with changes in ocean transport convergence on timescales longer than 10 years (e.g. Knight et al., 2005; Kim et al., 2018; Oelsmann et al., 2020; Fang et al., 2021). Furthermore, other studies have shown that the dynamical adjustment timescale of the Atlantic circulation in response to the volcanic forcing can be up to 15 years (e.g. Mignot et al., 2011; Swingedouw et al., 2015). Therefore, the predictions might be too short to reflect an impact due to the volcanic forcing on the AMV.

To further explore if the AMV response could be explained by changes in North Atlantic Ocean dynamics, we investigate the changes in the mixed layer depth and the AMOC (i.e. precursors to AMV changes). First we consider the multi-model and multi-eruption composites for the mixed layer depth in the subpolar North Atlantic. Figure 15a–d show a significant enhancement of the deep mixed layer in late winter and early spring (when the mixed layer depth in this region attains its maximum) in the first 3 years following the eruptions, with a maximum in the Labrador Sea (Fig. 15a–d). The mixed layer depth is usually interpreted as a proxy for deep convection, and in many models enhanced deep mixing in the subpolar North Atlantic is associated with strengthening of the AMOC a few years later (e.g. Dong and Sutton, 2005; Ortega et al., 2011). Indeed, Fig. 15e–h show that the mean AMOC stream function also experiences a small but significant strengthening in response to the volcanic forcing during years 2–9 after the eruption. The sign of the response is consistent for all the individual eruptions and the magnitude of both the mixed layer deepening and the AMOC strengthening seems to depend on the magnitude of the volcanic forcing, with responses in order of decreasing magnitude for Pinatubo, Agung and El Chichón.

Figure 15Multi-model and multi-eruption composites response for the mixed layer (February–March–April) depth (m) for years 1–3 and overturning stream function (Sv) years 2–9 to the volcanic eruptions. Stippling indicates statistically significant anomalies (see the Methods section).

We now explore the temporal evolution of the responses in the whole ensemble for two indices of the mixed layer depth and AMOC by averaging over the purple boxes in Fig. 15d and h, respectively. Figure 16 shows that the response of these indices has a large spread among the models, with HadGEM3-GC31-MM and CMCC-CM2-SR5 dominating the multi-model signal and the rest of the models not simulating a consistent or a significant response. Hence, Fig. 15 should be interpreted with care as the result is dominated by these models.

Figure 16Volcanic response (DCPP-A minus DCPP-C) of the mixed layer depth in the subpolar North Atlantic (50–65° N, 60–35° W) in February–March–April (top row) and the annual mean overturning stream function averaged over 30–50° N and 1000–2000 m (bottom row). Filled circles and triangles indicate statistically significant differences (see the Methods section). The ensemble mean for each model and the multi-model mean are shown. The shading is the multi-model member spread calculated as the 10th and 90th percentiles of the entire ensemble. The mixed layer depth for CESM1-1-CAM5-CMIP5 was not available and therefore could not be included.

Download

To further explore the origin of the inter-model differences in the AMOC response we look at the background density stratification conditions in each model to infer their preconditioning role in convection. Figure 17a–c show the change in the multi-eruption composite profiles of temperature, salinity and density in the Labrador Sea due to the effect of the volcanic eruptions in the first three winters (DJF) after their occurrence. We focus on DJF, which is a couple of months before the peak season for convection, because it includes both the preconditioning signal and the response to the vertical mixing, which hinders the interpretation of the results. CMCC-CM2-SR5, HadGEM3-GC31-MM and IPSL-CM6A-LR show significant changes in temperature, salinity and density in response to the volcanic forcing, which mostly imply a cooling, a salinification and a densification of the upper ocean levels. The simulated increase in upper density erodes the mean stratification and is thus consistent with a subsequent enhancement in convection. As previously shown, this is indeed the case for the CMCC-CM2-SR5 and HadGEM3-GC31-MM models (Fig. 16), for which the specific processes at play might differ as salinity changes seem to dominate the density changes in CMCC-CM2-SR5, while for HadGEM3-GC31-MM temperature changes seem to be more important. The reason for the lack of responsiveness of the mixed layer depth in IPSL-CM6A-LR, despite the significant surface densification exerted by the volcanic forcing, is the very strongly stratified mean state density conditions that it presents in the region, which are linked to very fresh upper-ocean conditions (Fig. 17d–f). In the rest of the models no significant changes in temperature or salinity are observed in response to the eruptions, which explains why neither the mixed layer depth nor the AMOC show a significant response.

Figure 17Multi-eruption mean response of the DJF Labrador Sea (50–65° N, 60–35° W) potential temperature, salinity and potential density profiles in years 1–3 after the eruptions (a–c) and absolute DJF Labrador Sea potential temperature, salinity and potential density (sigma1) profiles in years 1–3 after the eruptions in the predictions without volcanic forcing (d–f). Crosses indicate statistically significant differences (see the Methods section). The ensemble mean for each model and the multi-model mean are shown. The shading is the multi-model member spread calculated as the 10th and 90th percentiles.

Download