Fingerprint of volcanic forcing on the ENSO–Indian monsoon coupling

Phase coherence analysis

The ENSO and IM relationship is known to exhibit nonstationary behavior (9, 29), and its time-dependent evolution is in turn related to variations in both the amplitude and phase of the two oscillatory systems (30). In particular, phase coherence information is known to be critical for quantifying nonlinear interrelations between irregular, quasi-periodic, and oscillatory signals that cannot be inferred from the conventional linear correlation and wavelet approaches (30–33). We quantify the influence of LVEs on ENSO-IM coupling using phase synchronization between the monthly time series of phases of ENSO and IM computed from the PMIP3 last-millennium simulations (see Eq. 1, Materials and Methods; fig. S1). It can be seen that several epochs of ENSO-IM phase coherence, indicated by plateaus in the monthly time series of the phase difference between ENSO and IM, coincide with LVEs (Fig. 1A).

Event synchronization between phase-coherent epochs and LVEs

The coincidence between periods of ENSO-IM coupling and LVEs is also evidenced in observational records of ENSO and IM variations since the 19th century (fig. S2 and table S2). We applied event synchronization (ES) analysis [see Eq. 2, Materials and Methods; (34)] on the event series of ENSO-IM phase coherence and LVEs to quantify the coincidence of volcanic-induced ENSO-IM coupling. We uncover several epochs of ENSO-IM phase coherence. Our analysis suggests that LVEs tend to be followed by periods of ENSO-IM phase coherence, with a statistically significant ES strength (P < 0.02; 1000 event surrogates prepared by random shuffling) between the periods of phase-locking and LVEs.

Statistical significance of phase coherence analysis using recurrence-based twin surrogates

Owing to the complex, nonlinear and nonstationary connection between LVEs and ENSO-IM coupling, we performed a statistical significance test based on 5000 recurrence-based twin surrogates of ENSO and IM (see section S1), which is presented in fig. S2. The null distribution of the phase difference between ENSO and IM variations estimated from the 5000 twin surrogates is compared with the distribution of the actual phase difference based on the last-millennium simulations. We observe peaks at 0 and 2π in the distribution of ENSO-IM phase difference in contrast to the flat distribution from the null model, indicative of the signatures of LVE-induced ENSO-IM phase coherence (fig. S2, lower right panel). A similar distribution of ENSO-IM phase coherence analysis is also revealed in the instrumental records for the period 1871–2016 (fig. S2, lower left panel).

Bayesian analysis of LVE-induced ENSO-IM coupling

Compelling evidence for LVE-induced ENSO-IM coupling also comes from analyzing 25 paleoclimatic records (see datasets, table S1): 14 ENSO proxies and 11 IM proxies covering the last millennium (Fig. 1B). Conditional probability-based Bayesian analysis (35) of ENSO-IM coincidence in paleoclimatic reconstructions is performed, subject to (i) ENSO, (ii) ENSO and LVE in year −1, (iii) ENSO and Pacific Decadal Oscillation (PDO) (36), and (iv) ENSO and PDO and LVE in year −1 (see Eqs. 3 to 8). We find that 94 of 154 combinations showed that the coincidence of (El Niño and IM deficit)/(La Niña and IM excess) during a particular year was substantially higher when preceded by an LVE in the previous year, relative to the unconditioned probabilities (Fig. 1B). LVE-induced large-scale climatic impacts are also consistent with anomalous surface cooling over the Eurasian landmass evidenced from the paleoclimatic reconstructions of surface air temperature variations during the last millennium (Fig. 1C). Repeating the Bayesian analysis for different lead times after LVEs (Eq. 5) suggests that the probability of ENSO-IM coincidence decreases with each passing year until +4 years relative to LVE (fig. S3). Several epochs of high anti-correlation between ENSO and IM, e.g., 1401–1411, 1630–1640, 1715–1725, and 1888–1898, were coincident with high volcanic activity in paleoclimatic reconstructions [table S2; (12)]. These epochs can also be seen as plateaus in the ENSO-IM phase coherence analysis (Fig. 1A). In particular, the decade of 1630–1640 is reported to have experienced prolonged famines caused by successive crop failures following monsoonal breakdowns, leading to around 7.4 million deaths during 1630–1632 (37).

LVE-induced modulations of the angular frequency of ENSO and IM oscillations

Using the last-millennium simulations, we computed the time series of the angular frequency of ENSO and IM oscillations by taking the time derivative of their instantaneous phases (see Materials and Methods). The angular frequency of the ENSO was later segregated into two categories: (i) cases when ENSO events occurred in a 4-year window after an eruption (Volc) and (ii) cases when ENSO events were not preceded by LVEs during the preceding 4 years (NoVolc). The sampling distributions corresponding to the Volc and NoVolc categories were estimated using bootstrapping. A quantile-quantile (Q-Q) plot of the angular frequency of ENSO between the Volc and NoVolc cases shows an enhanced probability of higher angular frequency of ENSO for Volc events, as compared to NoVolc events (Fig. 2A). A Kolmogorov-Smirnov (K-S) test (35) indicates that the difference in the phase speed of ENSO between the Volc and NoVolc distributions is significant (P = 0.097) (see Materials and Methods), pointing to the role of LVEs in increasing the phase speed of ENSO. A similar Q-Q analysis was performed for the distributions of angular frequency of IM from the Volc and NoVolc categories (fig. S4) from the observational data during the 20th century. We note enhanced probability of lower angular frequency of the IM for Volc events as compared to NoVolc events. A K-S test for the observational data indicates that the difference in the phase speed of ENSO (P = 0.0001) and IM (P = 0.096) oscillations between the Volc and NoVolc distributions is statistically significant, thereby pointing to LVE-induced increases in angular frequency of ENSO and enhanced ENSO-IM phase coherence.

Role of background state on LVE-induced ENSO-IM coupling

While LVEs tend to increase the likelihood of triggering El Niño events (17), it is unclear whether every LVE would lead to an increase in the probability of ENSO-IM phase synchronization. A case in point is the famous Krakatoa eruption in 1883—one of the strongest volcanic events in the last two centuries, which was not followed by an El Niño in the subsequent year (38). While the timing of the Krakatoa eruption in August 1883 was possibly too late in the year to affect the equatorial Pacific Ocean in the first winter, it must be noted that the event also coincided with a slow decadal-scale transitioning of the Pacific SST anomalies into a La Niña–like cold phase, associated with the PDO, as evidenced from SST reconstructions based on a combination of coral and tree-ring proxy data together with instrumental records (39, 40). The decades after the Krakatoa 1883 eruption, i.e., the late 1880s to 1900s, coincide with one of the strongest epochs of ENSO-IM coupling during the 19th and 20th centuries (30). Comparing the magnitudes of ENSO-IM coincidence probabilities, we uncover that the probability of LVE-induced ENSO-IM coupling (Fig. 1B) is even more pronounced when the phase of ENSO is supported by the phase of the PDO (fig. S5, A and B), as they evolve in the background of the tropical Indo-Pacific coupled system. While the ENSO-IM coupling tends to be favored by the phases of ENSO and PDO coevolving in the same direction (left panels of fig. S5), we also note that the coincidence probabilities of ENSO-IM coupling are further enhanced until year +4 following an LVE (right panels of figs. S5 and S6).

Ocean initial conditions, strength of volcanic forcing, and ENSO evolution

We further evaluate the evolution of ENSO and its dependence on ocean initial conditions (ICs) and strength of volcanic radiative forcing (VRF), by conducting targeted ensemble climate sensitivity experiments using the IITM-ESM (28). The large internal variability of the climate system necessitates the use of a large number of ensembles for quantifying the forced response [fig. S7; (41, 42)]. Before proceeding to the discussion on the forced response, we shall first examine the role of ocean ICs on the transient evolution of Niño3 SST variations in the preindustrial control run (without volcanic forcing) of the IITM-ESM.

For the ensemble realizations, we first identified 100 oceanic states using the annual mean Niño3 SST anomalies that were randomly drawn from the 500-year preindustrial control run of the IITM-ESM (see Materials and Methods). We further grouped the 100 oceanic states into three categories of ICs based on the January SST anomalies, viz., 24 warm ICs (January Niño3 SST anomaly > 0.5°C), 29 cold ICs (January Niño3 SST anomaly < −0.5°C), and 47 neutral ICs (−0.5°C < January Niño3 SST anomaly < 0.5°C) (see Materials and Methods). All the targeted ensemble experiments start from January 1st, corresponding to the three categories of initial states of the ocean-atmosphere coupled system. An important aspect relevant to our study is to understand the transient evolution of ENSO with and without LVEs. To evaluate the transient evolution of ENSO in the absence of volcanic forcing, we first considered 3-year segments of monthly Niño3 SST anomalies of the preindustrial control run starting from the different ICs. Henceforth, we shall refer to the 3-year segments of monthly Niño3 SST anomalies starting from the 100 initial states of the preindustrial control run as the reference (REF) transient simulation, which serves as the baseline for comparing the transient evolution of ENSO in the volcanic forcing experiments.

To understand the transient evolution of ENSO anomalies in the absence of volcanic forcing, we first examined the means ± standard error of the monthly Niño3 SST anomalies of the REF ensembles from January of year −0 through December of year −2, for the three categories of ICs (Fig. 3). We note a general tendency for the Niño3 SST anomalies of REF starting from the warm (cold) ICs to transition into cooler-than-normal (warmer-than-normal) anomalies, respectively, after about 2 years (24 months), with considerable spread across the ensemble members. Likewise, the Niño3 SST anomalies of REF corresponding to the neutral IC ensemble tend to attain near-normal values after about 2 years. It is also interesting to note from Fig. 3 that the Niño3 SST anomalies corresponding to the warm and neutral ICs of the REF experiment tend to continue in warmer-than-normal states at the end of the first year (12 months), whereas the Niño3 SST anomalies corresponding to the cold ICs tend to continue in cooler-than-normal states at the end of the first year.

For evaluating the behavior of ENSO-IM coupling under varying magnitudes of VRF, we performed a 100-member ensemble experiment (VRF1x) of the IITM-ESM forced by the Krakatoa VRF from year −0 (January 1883) through year −2 (December 1885). The vertically integrated aerosol optical depth (AOD) at 550 nm during the Krakatoa eruption (Fig. 4A) is a good measure of the magnitude of the global VRF. The VRF1x ensemble (100-member) simulations were started from the same set of January 1st ICs corresponding to the three categories (cold, neutral, and warm) of initial states, and the model was integrated for 3 years (see Materials and Methods).

To assess the LVE-induced modulations of ENSO, we examined the means ± standard error of the ΔNiño3 SST variations (i.e., the member-to-member difference of Niño3 SST between the VRF1x and REF experiments) from January 1883 through December 1885, starting from the neutral, warm, and cold ICs, respectively (Fig. 4). Since the VRF1x and REF ensemble members start from identical ICs, the member-to-member difference (VRF1x − REF) provides a measure of the Niño3 SST response to the external volcanic forcing. While it is seen that the ΔNiño3 variations for the neutral IC ensembles exhibit relatively smaller fluctuations around 0°C, it is important to note that the ΔNiño3 variations for the warm IC ensembles first rapidly transition to negative values (mean ~ −1.08°C) during the winter (December-January-February) of 1883–1884 and then subsequently undergo a transition to positive values (mean ~ 0.3°C) during the winter of 1884–1885. For the cold IC ensembles, we see an opposite behavior with the ΔNiño3 variations first transitioning rapidly to positive values (mean ~ 0.5°C) around the winter of 1883–1884 and later attaining near-neutral values around the winter of 1884–1885. The rapid transitioning of the cold IC ensembles to positive anomalies of ΔNiño3 during the first winter (i.e., 1883–1884), as compared to the warm ICs and neutral ICs, seems to agree with the findings of (18) for high-latitude eruptions. The fast changes in the phase of ΔNiño3 seen in Fig. 4 are consistent with an increase in the angular frequency of ENSO in response to external forcing from LVEs, as discussed earlier in Fig. 2.

Of the 100 members, we find that 41 members transitioned to warmer-than-normal (ΔNiño3 > 0.5°C), 25 members transitioned to cooler-than-normal (ΔNiño3 < −0.5°C), and the remaining 34 members transitioned to near-neutral (−0.5°C < ΔNiño3 < 0.5°C) states, during the winter of 1884–1885, thus indicating a significantly higher probability of transitioning toward an El Niño in the year following the eruption (fig. S7A). We further note that of the 41 warmer-than-normal transitions, 13 members started from cold ICs, 20 members started from neutral ICs, and 8 members started from warm ICs, which suggests that cold ICs tend to have a higher probability (0.45 = 13/29) of transitioning to warmer-than-normal Niño3 anomalies in the year following the eruption, as compared to the probability (0.33 = 8/24) of warm ICs, which is in agreement with the findings of (18). Note that the 100 ICs consist of 29 cold ICs, 47 neutral ICs, and 24 warm ICs (fig. S7A).

In addition, we performed two more sets of short (3 years) ensemble simulations by varying the magnitude of the Krakatoa VRF. In the first set (VRF4x), we quadrupled the magnitude of volcanic forcing, and for the second set (VRF0.25x), the magnitude was reduced to one-fourth. It must also be mentioned, here, that we use different sets of ICs for the VRF4x and VRF0.25x ensemble simulations. The VRF4x experiment is a 25-member ensemble simulation starting from the same ICs as those of the 25 cases that transitioned to cooler-than-normal anomalies (i.e., Niño3 anomaly < −0.5°C) during the winter of 1884–1885 in the VRF1x simulation (fig. S7B). Note that the abovementioned 25 realizations of VRF1x can be traced back to the January 1st ICs corresponding to the 7 cold ICs, 10 neutral ICs, and 8 warm ICs, respectively (fig. S7, A and B). A similar approach is used for the VRF0.25x experiment, except for a 41-member ensemble simulation in this case. The 41-member ensembles in the VRF0.25x start from the same ICs as those of the 41 cases that transitioned to warmer-than-normal anomalies (i.e., Niño3 anomaly >0.5°C) during the winter of 1884–1885 in the VRF1x simulation (fig. S7C). Note that the abovementioned 41 realizations of VRF1x can be traced back to the January 1st ICs corresponding to the 13 cold ICs, 20 neutral ICs, and 8 warm ICs (fig. S7, A and C).

To examine the effect of quadrupled volcanic forcing on the transient ENSO response, we compared the time evolution of ΔNiño3 SST variations between the 25 ensemble members of VRF4x and VRF1x experiments (Fig. 5A). Note that the 25 members of VRF4x are compared with those realizations of VRF1x that started from the same 25 ICs (fig. S7B). The member-to-member comparison between the ensembles of VRF4x and VRF1x allows us to distinguish the ΔNiño3 SST response to changes in the magnitude of volcanic forcing from the internal dynamics of the coupled system. We find that the ΔNiño3 variations for the VRF4x ensembles exhibit a rapid transition to negative values from August 1883 through July 1884 and later evolve to positive anomalies from August 1884 through December 1884. Anomalous cooling (mean ~ −1.04°C) of the ΔNiño3 response in the VRF4x 25-member ensembles is prominently seen around December 1883, while the warm anomaly (mean ~ 0.7°C) peaks around November 1884 (Fig. 5A). In contrast, the VRF1x ensembles show an anomalous cooling in the ΔNiño3 response from September 1884 through August 1885, with a substantial cooling (mean ~ −0.83°C) around January 1885 (Fig. 5A). We further note that 5 of the 25 realizations in VRF4x transitioned to warmer-than-normal conditions (Niño3 anomaly > +0.5°C), 9 members transitioned to neutral conditions (−0.5°C < Niño3 anomaly < +0.5°C), and 11 members transitioned to cooler-than-normal conditions (Niño3 anomaly < −0.5°C) during the winter of 1884–1885 (fig. S7B). The above results suggest that extremely strong volcanic forcing (e.g., quadrupled Krakatoa) can significantly enhance not only the angular frequency of ENSO but also the probability of occurrence of El Niño during the latter half of the second year, following the eruption.

Furthermore, we also compared the ΔNiño3 SST variations between the 41 ensemble members of VRF0.25x and VRF1x, which start from the same ICs (fig. S7C), to understand the transient ENSO response to a reduced volcanic forcing (Fig. 5B). It can be seen that the ΔNiño3 variations in the 41-member VRF1x ensembles transition to positive anomalies from July 1884 through May 1885, with a pronounced warming (mean ~ 1.29°C) around December 1884 (Fig. 5B). On the other hand, the anomalous warming of the ΔNiño3 response is muted in the VRF0.25x ensembles during the latter half of second year and beyond (Fig. 5B), which suggests that decreasing the strength of the volcanic forcing tends to weaken the simulated El Niño–like response by the end of the second year.

Physical insights into the effects of LVEs on the response of the tropical atmosphere-ocean coupled system can be deduced by examining the spatial composite maps of the differences in surface temperature, depth of the 20°C isotherm (D20), precipitation, and net surface shortwave radiation fluxes between the VRF1x and REF ensembles (Fig. 6). The composites are based on the average of the member-to-member difference between the 100-member ensembles of the VRF1x and REF realizations. The ΔSST response to LVE in the year following the eruption shows an El Niño–like pattern of anomalous warming in the equatorial eastern-central Pacific and cooling in the west Pacific (Fig. 6A). In addition, anomalous cooling can be noticed over the Northern Hemispheric (NH) land region in the surface air temperature response (Fig. 6A). The ΔD20 anomalies (Fig. 6B) show anomalous deepening of the thermocline in the equatorial eastern Pacific, indicating an El Niño–like response in the year following an eruption. It is also interesting to note an anomalous deepening of thermocline in the equatorial eastern Indian Ocean (Fig. 6B), which is reminiscent of a negative Indian Ocean Dipole (IOD) (see Supplementary Materials and Methods, section S3). While the thermocline deepening is prominent in the eastern equatorial Indian Ocean (Fig. 6B), the overlying SST anomalies are relatively weaker in magnitude (Fig. 6C) probably because of the volcanic cooling. However, the negative IOD-like pattern emerges more clearly in the relative SST anomalies after subtracting the mean value of the tropical Indo-Pacific SST (Fig. 6E).

The precipitation response in the year following the eruption shows decrease of the boreal summer monsoon precipitation over South and Southeast Asia and over the tropical west Pacific to the north of equator (Fig. 6C). An increase of precipitation is noticed south of the equator over the Pacific Ocean and the Indian Ocean (Fig. 6C). The pattern of a dry (wet) precipitation anomaly to the north (south) of the equator is also seen over the West African monsoon region, eastern Pacific, and adjoining areas. The imposition of volcanic forcing leads to a pronounced decrease of the net shortwave radiative fluxes at the surface by nearly 5 to 10 W m⁻² over the NH, with smaller reductions over the Southern Hemisphere, thus giving rise to an inter-hemispheric contrast in the net shortwave radiative fluxes (Fig. 6D). Earlier studies have pointed to the effects of differential cooling of the NH in triggering a southward shift of the inter-tropical convergence zone (ITCZ) and creating an El Niño–like response (24–26), following volcanic eruptions. On the basis of our overall understanding, as well as the results shown in Fig. 6, it is suggested that the volcanic-induced cooling over the NH extra-tropics induces a southward shift of the ITCZ and causes an El Niño–like response in the tropical Pacific, as well as a negative IOD-like response in the tropical Indian Ocean (Fig. 6E, Supplementary Materials and Methods, section S3). An increase of precipitation south of the equator, resulting from the anomalous southward shift of the ITCZ and the negative IOD-like response (Fig. 6, B and E), tends to suppress monsoon precipitation over South and Southeast Asia, through a weakening of the regional summer monsoon Hadley-type overturning circulations (see Supplementary Materials and Methods, section S3).

Prompted by the contrasting ΔNiño3 SST responses between the VRF4x and VRF1x ensembles (Fig. 5A), we additionally examined the composite spatial maps of the member-to-member difference (relative to REF) in surface temperature and precipitation for both the VRF4x and VRF1x ensembles. The composited maps of the member-to-member difference in surface temperature for (VRF4x – REF) and (VRF1x – REF), based on the (Mem#25) ensembles, are shown in fig. S8 (A and B, respectively). It is to be noted that all the members of the (Mem#25) ensembles of VRF4x, VRF1x, and REF (fig. S8) originate from the same ICs. Further, it can be seen that the SST response for (VRF1x – REF) in fig. S8B shows cold anomalies in the central-eastern Pacific, which is consistent with the transitioning of the (Mem#25) ensembles of VRF1x to cooler-than-normal (ΔNiño3 < −0.5°C) anomalies during the winter of 1884–1885 (fig. S7A). In contrast, the SST response for (VRF4x – REF) in fig. S8A shows an El Niño–like pattern with anomalous warming in the equatorial eastern Pacific and cooling in the western Pacific. It is also seen that the anomalous cooling is stronger in (VRF4x – REF) as compared to (VRF1x – REF), over most of the NH extra-tropical land and oceanic regions. We also reckon that the warm anomaly over the region of northeast Asia (fig. S8A) may be part of the natural internal variability over the mid and high latitudes and that the relatively weaker NH cooling in the (Mem#25) ensembles, as compared to the 100 members, could explain the different ENSO response (18).

The composited maps of the member-to-member difference in precipitation for (VRF4x – REF) and (VRF1x – REF), based on the (Mem#25) ensembles, are shown in fig. S8 (C and D, respectively). We note a strong suppression of precipitation over the summer monsoon regions of South and Southeast Asia, tropical west Pacific to the north of the equator in (VRF4x − REF), while the increased precipitation to the south of the equator over the central-eastern Pacific, Indian, and Atlantic oceans is associated with the southward shift of the ITCZ (fig. S8C). Also seen in fig. S8C is the strong decrease of summer precipitation to the north of equator over the eastern Pacific, Atlantic, and the African monsoon region. In contrast, the precipitation response in (VRF1x − REF) exhibits above-normal precipitation over Southeast Asia (fig. S8D), which is consistent with the corresponding SST response shown in fig. S8B.

In addition, we examined the effects of reduced volcanic forcing on the tropical ocean-atmosphere system by comparing the (VRF0.25x – REF) and (VRF1x – REF), based on the (Mem#41) ensembles (fig. S8, E to H). We note from fig. S8 (E and F) that the anomalous SST response (warm in the equatorial central-eastern Pacific and cold in the tropical western Pacific), as well as the anomalous cooling over the NH extra-tropics are considerably damped in (VRF0.25x – REF) as compared to (VRF1x – REF). Associated with the damped surface temperature response, a slight increase in monsoon precipitation over South and Southeast Asia is seen in (VRF0.25x – REF) as compared to (VRF1x – REF) (fig. S8, G and H). In short, the findings from the large-ensemble targeted sensitivity experiments presented above suggest that the global and tropical ocean-atmosphere coupled response to LVEs serves as an effective mechanism for coupling the ENSO and IM oscillatory systems.