Research ArticleGEOLOGY

Fingerprint of volcanic forcing on the ENSO–Indian monsoon coupling

See allHide authors and affiliations

Science Advances  18 Sep 2020:
Vol. 6, no. 38, eaba8164
DOI: 10.1126/sciadv.aba8164

Abstract

Coupling of the El Niño–Southern Oscillation (ENSO) and Indian monsoon (IM) is central to seasonal summer monsoon rainfall predictions over the Indian subcontinent, although a nonstationary relationship between the two nonlinear phenomena can limit seasonal predictability. Radiative effects of volcanic aerosols injected into the stratosphere during large volcanic eruptions (LVEs) tend to alter ENSO evolution; however, their impact on ENSO-IM coupling remains unclear. Here, we investigate how LVEs influence the nonlinear behavior of the ENSO and IM dynamical systems using historical data, 25 paleoclimate reconstructions, last-millennium climate simulations, large-ensemble targeted climate sensitivity experiments, and advanced analysis techniques. Our findings show that LVEs promote a significantly enhanced phase-synchronization of the ENSO and IM oscillations, due to an increase in the angular frequency of ENSO. The results also shed innovative insights into the physical mechanism underlying the LVE-induced enhancement of ENSO-IM coupling and strengthen the prospects for improved seasonal monsoon predictions.

INTRODUCTION

With nearly 50% of the net cultivated area of India being rain-fed (1), accurate forecasts of the Indian monsoon (IM) are of paramount importance for the agrarian economy of the country (2). Year-to-year variations in the IM are closely tied to large-scale patterns of sea surface temperature (SST) and sea level pressure variations in the tropical Pacific Ocean associated with the El Niño–Southern Oscillation (ENSO) phenomenon, wherein periods of anomalously dry IM and wet IM coincide with El Niño and La Niña conditions, respectively (310).

Large volcanic eruptions (LVEs) have a strong impact on Earth’s climate variability from interannual to decadal time scales (11, 12). Volcanic ash and gaseous matter from LVEs injected into the stratosphere alter the global mean surface temperature through backscatter and absorption of shortwave radiation, leading to changes in atmospheric and ocean circulation and the global hydrological cycle (13, 14). Abrupt changes in global climatic patterns induced by LVEs can have serious societal implications, e.g., the Mt. Tambora eruption in Sumbawa, Indonesia, in 1815, resulted in the “year without a summer” and caused a severe frost-induced socioeconomic breakdown in much of Europe and North America (15). Within the broader climatic context, recent studies have focused on the modulation of the ENSO by LVEs (1620), as well as the global and Asian monsoon precipitation response to LVEs (2124). Although we now have a better understanding of the impact of LVEs on the ENSO system and their ability to trigger El Niños, in particular (19, 2426), there has been little or no focus on the impact of LVEs on the coupling between ENSO and IM. Here, we investigate how LVEs, besides directly influencing the ENSO and IM, might influence the interrelationship between the two systems. Our analysis is based on historical observational data of SST variations in the tropical Pacific Niño3 region, IM rainfall, together with paleoclimate reconstructions of ENSO and IM variations during the last millennium based on multiple proxies (see datasets, table S1), last-millennium climate simulations from the Paleoclimate Modeling Intercomparison Project phase 3 (PMIP3) (27), and large-ensemble targeted climate model sensitivity experiments using the Indian Institute of Tropical Meteorology Earth System Model (IITM-ESM) (28).

RESULTS

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 (3033). 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).

Fig. 1 LVE-induced ENSO-IM phase coherence during the last millennium.

(A) Phase coherence analysis of ENSO and IM using Institut Pierre Simon Laplace (IPSL) PMIP3 last-millennium simulations (850–1850 AD). Time series of phase difference between the instantaneous phases of ENSO and IM (top). The green shading shows periods of phase coupling between ENSO and IM. Time series of VRF (W m−2) during the last millennium based on (12) (bottom). (B) Bayesian analysis—probability of ENSO-IM phase coupling for year (0) conditioned on ENSO (0) (left) and ENSO (0) and Volcano (−1) (right). The paleoclimate reconstructions of 14 ENSO indices (4659) and 11 IM indices are based on multiple datasets (see table S1, EN1,…EN14/ IM1,…IM11). (C) Composite of Northern Hemispheric (NH) surface air temperature (SAT) anomalies based on LVEs during the last millennium. The composite is based on 38 cases of LVEs (12, 17). The SAT anomalies are based on the tree-ring reconstructions (60). The composite demonstrates prominent NH cooling following LVEs.

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.

Fig. 2 LVE-induced modulations of the angular frequency of ENSO and IM oscillations during the last millennium (850–1850).

Q-Q plot showing angular frequency of (A) ENSO and (B) IM oscillations based on the volcanic and no volcanic distributions. The volcanic distributions are bootstrapped from year (0) to year (4) relative to volcanic eruption. Angular frequencies are obtained as the first derivative of the instantaneous phase computed using phase coherence analysis of monthly indices of Niño3 and IM precipitation from the IPSL PMIP3 last-millennium simulation. The oscillations of ENSO evolve slowly relative to the IM oscillations, so that increases in the angular frequency of ENSO tend to enhance 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.

Fig. 3 Transient evolution of ENSO in the REF ensembles (without volcanic forcing) starting from different oceanic initial states.

Temporal evolution of the variations of Niño3 SST anomalies from the 100-member REF ensembles without volcanic forcing (i.e., based on the preindustrial control simulation) shown from January of year −0 through December of year −2. The monthly Niño3 SST anomalies of the REF ensembles are relative to the long-term baseline climatology of the preindustrial control simulation. The variations of the Niño3 SST anomalies are expressed as means ± standard error of the ensemble members. The 100 members consist of realizations starting from 47 neutral ICs, 24 warm ICs, and 29 cold ICs, respectively.

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).

Fig. 4 LVE-induced modulations of ENSO under different ICs.

Temporal evolution from January 1883 to December 1885 (A) Global mean stratospheric AOD. (B) LVE-induced ENSO anomaly expressed as means ± standard error of the ΔNiño3 variations, computed from the member-to-member difference between the VRF1x and REF ensembles. The top, middle, and bottom panels correspond to ensemble realizations starting from 47 neutral ICs, 24 warm ICs, and 29 cold ICs, respectively.

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.

Fig. 5 Dependence of ENSO response on the strength of LVE.

Temporal evolution of ΔNiño3 variations, expressed as means ± standard error, shown from January 1883 to December 1885. (A) The two panels of ΔNiño3 variations (orange color) are computed from the member-to-member difference of the 25-member cluster (Mem#25) corresponding to (VRF4x − REF) and (VRF1x − REF), respectively. (B) Same as (A), except for the 41-member cluster (Mem#41) corresponding to (VRF0.25x − REF) and (VRF1x − REF), respectively (gray color). The blue window indicates the period from June 1884 to February 1885.

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).

Fig. 6 Tropical Indo-Pacific and monsoon response to Krakatoa VRF.

Spatial maps of the composited difference (VRF1x − REF) constructed from the 100-member ensembles of the VRF1x and REF realizations. (A) SST + SAT over land (°C) anomalies. (B) Depth of the 20°C (D20) isotherm (m) anomalies. (C) Precipitation (mm day−1) anomalies. (D) Clear sky net shortwave radiative fluxes at the surface (W m−2) anomalies and (E) SST anomalies in the Indian Ocean relative to the mean value in the tropical Indo-Pacific region (40°E to120°E; 8°S to 8°N). The composites are based on the average of the difference between the 100-member ensembles of the VRF1x and REF realizations. Before averaging, we first computed the member-to-member difference between the VRF1x and REF ensembles starting from the same ICs. The surface temperature anomalies are averaged from June 1884 through February 1885 (blue window in Fig. 5), while the precipitation and surface shortwave radiation fluxes are averaged for the June to September 1884 boreal summer monsoon season. RSST, relative SST.

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−2 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 (2426), 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.

DISCUSSION

Nonstationarity in the relationship between ENSO and IM is recognized as a major bottleneck for advancing seasonal predictions of the IM (9, 43). The investigation by Maraun and Kurths (30) on the phase dynamics of ENSO and IM oscillatory systems revealed the occurrence of distinct epochs of phase coherence of ENSO and IM oscillations in the historical data during 1871–2003, which were unlikely to be a result of stochastic fluctuations. They further conjectured that the radiative effects of volcanic eruptions could potentially induce epochs of ENSO-IM phase coherence (30). However, fingerprinting the role of LVEs on ENSO-IM coupling is a complex problem given the large internal variability of the climate system. These considerations inspired us to investigate the signatures of LVEs on the ENSO and IM systems using a panoply of long-term climate datasets, including the last-millennium climate simulations and multiple paleoclimate proxy reconstructions, and by using a variety of advanced techniques for analyzing phase coherence, ES, and probabilities of LVE-induced ENSO-IM coupling. We also performed a large ensemble of targeted climate sensitivity experiments using the IITM-ESM (28) to assess the role of the background state of the Pacific Ocean in determining the extent to which LVEs can influence ENSO-IM coupling.

Our results provide compelling evidence that shows that the radiative forcing by LVEs tends to significantly enhance the phase coherence of ENSO and IM oscillations, arising as a consequence of an increase in the angular frequency of ENSO in response to volcanic forcing. While the probability of LVE-induced ENSO-IM coupling is substantially enhanced from the year +1 through year +4 following an LVE, it is noted that extremely strong volcanic eruptions (e.g., quadrupled Krakatoa) significantly enhance the probability of occurrence of El Niño–like conditions in the year following the eruption. In addition, the probability of coincidence of {El Niño and IM drought} or {La Niña and IM excess} following an LVE is higher when the phases of ENSO and PDO are in the same direction. Results from the IITM-ESM large-ensemble climate sensitivity experiments indicate the occurrence of widespread decrease in summer monsoon rainfall over South and Southeast Asia following LVEs. It is noted from the model simulations that the suppression of monsoon rainfall over South and Southeast Asia following LVEs arises from the combined effects of an El Niño–like response in the tropical Pacific, strong cooling over the NH land regions, negative IOD-like conditions in the tropical Indian Ocean, and a southward shift of the ITCZ (Fig. 6; Supplementary Materials and Methods, section S3).

That volcanic eruptions could potentially trigger El Niño–like conditions had been widely investigated [Supplementary Materials and Methods, section S4; (1620, 2426)], but whether they could also force coupling of ENSO and IM had remained unclear, so far. Our result that LVEs lead to increased synchronization of the ENSO and the IM has an added significance as it points to an additional source of predictability of seasonal monsoon rainfall in the years following a volcanic event. Furthermore, the present findings on how the ENSO-IM coupling is modulated via volcanism can help in improving the state-of-the-art climate models and particularly modeling the IM more reliably in a changing climate and even more so for assessing the regional implications of geo-engineering experiments.

MATERIALS AND METHODS

Datasets

The observed and reconstructed datasets used in our analysis are summarized in table S1. These include gridded historical observational data of SST (44) and rainfall over India (45), as well as paleoclimate reconstructions of 14 ENSO (4659) and 11 IM indices (table S1) based on multiple proxies. The datasets are quoted in the manuscript as (table S1, Index). For example (table S1, HIM) refers to the observed IM rainfall dataset for the historical period 1871–2016. Information about the paleoclimate reconstructions of 14 ENSO (4659), 4 PDO (table S1, PDO1,…PDO4), and 11 IM proxies (table S1, IM1, IM2,…IM11) during the last millennium is tabulated in table S1. We also analyzed reconstructed paleoclimatic surface air temperature anomalies for the last millennium (60). For the large-ensemble targeted climate sensitivity experiments with the IITM-ESM (28), we used the stratospheric AOD forcing from CMIP6 (Supplementary Materials and Methods, section S2). In addition, we have analyzed the last-millennium climate simulations from the PMIP3 (27).

Phase coherence

Phase coherence, which was first observed by C. Huygens in 1665, is used to study the synchrony of periodic signals, i.e., when the difference in phases of two oscillators is approximately constant. More recently, phase coherence was shown to occur in high-dimensional chaotic signals as well, where a pair of chaotic oscillators show strongly correlated phases with almost entirely uncorrelated amplitudes (32). Quantitatively, phase coherence depends on the phase φ(t) of a self-sustained oscillatory signal x(t) that encodes the relative location of the signal within its oscillation cycle such that each cycle φ(t) increases monotonically from 0 to 2π. Phase coherence analysis documented by Maraun and Kurths (30) showed that post-1980 decoupling of the ENSO and IM (9) was symptomatic of a “2:1 phase relation”; i.e., an IM completes a cycle in half the time of an ENSO cycle, which is also observable in earlier time periods, viz., 1908–1921 and 1935–1943. Without necessarily assuming chaoticity, this powerful extension to aperiodic oscillations enables us to use phase coherence to quantify interrelations between irregular, quasi-periodic, and oscillatory signals, as is typically seen in climatic indices. There are clear advantages in doing so: (i) the 2:1 phase relation reported by Maraun and Kurths (30) is not observable with linear cross-correlation or wavelet analysis; (ii) instead of a coarse-grained sliding window–based time evolution (9), we obtain an instantaneous relation at the original time scale of the signal; and (iii) since the phase is a universal property of all oscillators, it can be used to study the coupling between structurally different physical systems. A recent study exhibited the use of phase coherence between the precipitation dipole of South America and Southern Hemispheric circulation patterns and showed that monsoonal variability in the region is primarily controlled by the Rossby wave dynamics (32).

Phase coherence analysis

Phase coherence analysis (30) is used to compute the instantaneous phases of self-sustained oscillatory systems such as ENSO and IM. This analysis is performed using monthly time series of observed historical data of Niño3 (44) and IM (45) for the period 1871–2016 (see table S1, HEN and HIM), and also using the PMIP3 last-millennium simulation of the Institut Pierre Simon Laplace (IPSL) model (27). The time series of ENSO (i.e., SST anomaly averaged in the Niño3 region 5°N to 5°S, 150°W to 90°W) and IM (i.e., precipitation averaged over the IM core zone 74.5°E to 86.5°E and 16.5°N to 26.5°N) are averaged for each month from the last-millennium PMIP3 simulation of the IPSL climate model. Following (30), a low-pass Butterworth filter (cutoff frequency = 0.75 cycles per year, order 1) is applied to the monthly time series, so that intra-annual oscillations are filtered out from the time series of ENSO and IM (fig. S1). A second-order differencing of the time series is further applied followed by a Hilbert transform (Anal), which yields the real and imaginary components of the time series at every time step. The phase φ is computed from the real and imaginary components of the time series as followsϕ=tan1Imag (Anal)Real (Anal)(1)where Anal is the Hilbert transform of the time series.

Upon obtaining absolute instantaneous phases, the phase difference between Niño3 and IM is examined. The plateaus in the time series of the phase difference show regions of phase locking between ENSO and IM, which are indicated by green shading in Fig. 1A. The absolute instantaneous phases of Niño3 and IM are then used to compute the first difference of the respective time series as the angular frequency of ENSO and IM systems (Fig. 2 and fig. S4).

Identification of phase-locked periods

Identification of phase-locked periods is based on the distribution of ENSO-IM phase difference time series. The distribution shows a peak around zero. The selection criteria for the points are such that 12.5% of the periods captured on either side of the peak gives a good estimate of the plateaus in the phase difference plot. In addition, points on the trajectory in the complex plane that are near the center (radius around the center <0.00007 for ENSO and <0.00003 for IM) are excluded to avoid spurious oscillations of ENSO and IM systems in the identification. The extreme points are also excluded as they are not a part of the oscillatory system around a single attractor.

Event synchronization

ES is a technique to estimate synchronization between two event series (34). ES is carried out on the event series of volcanic eruptions (magnitude of global radiative forcing exceeding 0.25 W m−2) and epochs of phase coherence. The event series of phase coherence is coarse-grained to yearly resolution. We classified a year to be phase-coherent for the ENSO-IM system if at least 3 months of the year exhibit phase coherence. Using 100 randomly generated surrogates for the null model, a P value < 0.02 is obtained from the ES analysis. Following (34), the ES index Qτ for two simultaneously measured discrete univariate time series xn and yn, n = 1, …, N, with event times tix and tjy (i = 1, …, mx; j = 1, …, my, where mx and my are the number of events in xn and yn, respectively), is found asQτ=cτ(yx)+cτ(xy)(mx*my)(2)wherecτ(xy)=Σi=1mx  Σj=1my  Jijτ1if 0<tixtjyτijJijτ=1/2if tix=tjy0else,andτij=min{ti+1xtix,tixti1x,tj+1ytjy,tjytj1y}/2

Bayesian analysis

Bayesian analysis is used to diagnose the occurrence probability of El Niño–IM droughts and La Niña–IM excess, conditional to the occurrence of volcanic eruptions. We use paleoclimate reconstructions of 14 ENSO and 11 IM indices (see table S1, EN1,…EN14/IM1,…IM11) and forcing of volcanic eruptions (12) at an annual temporal resolution for the last millennium to evaluate the probability of ENSO-IM cooccurrence by conditioning with and without volcanic eruptions. Volcanic eruptions considered for this analysis are those events for which the magnitude of global radiative forcing exceeds 0.25 W m−2. The probability matrix of ENSO-IM coupling conditioned on (i) ENSO only and (ii) ENSO with volcanic eruption occurrence in the preceding year using the Bayesian analysis is shown in Fig. 1B.

To understand the role of additional factors influencing ENSO-IM cooccurrence in relation to LVEs, the probability of ENSO-IM coupling is diagnosed in such a way that volcano occurrence noted in the previous year is also conditioned on the background state of the tropical Indo-Pacific climate system. For this purpose, indices of PDO based on four paleoclimate reconstructions (see table S1, PDO1, PDO2, PDO3, and PDO4) are used for assessing the background state of the Indo-Pacific climate system. The matrices of the conditional probabilities of ENSO-IM cooccurrence conditioned on the state of PDO with and without LVEs are shown in Fig. 3A and figs. S5 and S6.

Conditional probability of ENSO-IM coupling

Conditional probability of ENSO-IM coupling is calculated by first extracting the common periods from the paleoclimatic time series of ENSO and IM. If there is no common period between the two ENSO and IM series, we get a null value of cooccurrence probability for that particular combination. The corresponding pixels are shown as white in the Bayesian analysis (Fig. 1B and figs. S3, S5, and S6). For the calculation of ENSO-IM cooccurrence probability, a count of the total number of droughts (monsoon excess) with and without El Niño (La Niña) is estimated. The conditional probability is then estimated asP(ENSOIM cooccurrenceENSO)=X+YZ+W(3)where X is the total count of drought events cooccurring with El Niño, Y is the total count of monsoon excess events cooccurring with La Niña, Z is the total count of El Niños, and W is the total count of La Niñas.

To check the volcanic link to ENSO-IM coupling, another condition is applied to the conditional probability, i.e., whether a volcanic eruption (magnitude >0.25 W m−2) occurred in the year before (Fig. 1B). The ENSO-IM cooccurrence probability conditioned on LVEs is then estimated asP(ENSOIM cooccurrence(ENSOVolcano (year=1 relative to cooccurrence)))=X+YZ+W(4)where X is the total count of drought events cooccurring with El Niño and eruption in the previous year, Y is the total count of monsoon excess events cooccurring with La Niña and eruption in the previous year, Z is the total count of El Niños with eruption in the previous year, and W is the total count of La Niñas with eruption in the previous year.

We also checked whether the eruption occurred at lags of 2, 3, and 4 years (fig. S3). The ENSO-IM cooccurrence probability conditioned on LVE at a lag of n years is estimated as follows:P(ENSOIM cooccurrence(ENSOVolcano (year=n relative to cooccurrence)))=X+YZ+W(5)where X is the total count of drought events cooccurring with El Niño and eruption n years before El Niño-drought events, Y is the total count of monsoon excess events cooccurring with La Niña and eruption n years before La Niña–monsoon excess events, Z is the total count of El Niños with eruption n years before El Niño, and W is the total count of La Niñas with eruption n years before La Niña.

It is noted that the probability of coupling declines as we move away from the volcanic eruption. The selection of events is followed in such a way that the time series of ENSO and IM are reproduced from raw data and also consistent with the source of these datasets. These time series are then normalized and converted into event series so that their normalized versions exceed ±0.5σ, for both ENSO and IM.

The conditioning on PDO also follows a similar way as adopted for ENSO (Fig. 2, A and B), except for the condition that the PDO warm (cold) phase coincides with El Niño (La Niña). It is estimated asP(ENSOIM cooccurrence(ENSOPDO))=X+YZ+W(6)where X is the total count of drought events cooccurring with El Niño and warm PDO, Y is the total count of monsoon excess events cooccurring with La Niña and cold PDO, Z is the total count of El Niños and warm PDO, and W is the total count of La Niñas and cold PDO.

The conditional ENSO-IM conditional probability conditioned on PDO and LVEs (fig. S5, A and B) is then estimated asP(ENSOIM cooccurrence(ENSOVolcano(year=1 relative to cooccurrence)PDO))=X+YZ+W(7)where X is the total count of drought events cooccurring with El Niño and warm PDO with eruption 1 year before El Niño–drought-warm PDO events, Y is the total count of monsoon excess events cooccurring with La Niña and cold PDO with eruption 1 year before La Niña–monsoon excess-cold PDO events, Z is the total count of El Niños and warm PDO with eruption 1 year before El Niño, and W is the total count of La Niñas and cold PDO with eruption 1 year before La Niña.

The probability of ENSO-IM cooccurrence conditioned on PDO and LVE is evaluated at 1-, 2-, 3-, and 4-year lags from the LVE (figs. S5 to S7). It is estimated asP(ENSOIM cooccurrence(ENSOVolcano(year=n relative to cooccurrence)PDO))=X+YZ+W(8)where X is the total count of drought events cooccurring with El Niño and warm PDO with eruption n years before El Niño–IM deficit-warm PDO events, Y is the total count of monsoon excess events cooccurring with La Niña and cold PDO with eruption n years before La Niña–monsoon excess-cold PDO events, Z is the total count of El Niños and warm PDO with eruption n years before El Niño, and W is the total count of La Niñas and cold PDO with eruption n years before La Niña.

Sensitivity of ENSO-IM coupling to the background state

Sensitivity of ENSO-IM coupling to the background state is evaluated by performing targeted ensemble simulations using the IITM-ESM. A brief description of the ICs for the targeted ensemble simulations and the volcanic forcing experiments is provided below. For the ensemble simulations, we have used 100 oceanic initial states that are identified from the 500-year preindustrial (PI-CTL) run of the IITM-ESM (28) and grouped them into three different categories (viz., warm IC, cold IC, and neutral IC) in a stepwise process, as described below.

In the first step, we use the annual mean SST anomalies and select 50 warm cases (annual mean Niño3 SST anomaly > 0.5°C) and 50 cold cases (annual mean Niño3 SST anomaly < −0.5°C). Here, it must be mentioned that the monthly and annual SST anomalies are relative to their respective baseline climatology constructed from the last 300 years of PI-CTL run. In the second step, we consider the January Niño3 SST anomalies for the 100 identified cases and group them into three categories, viz., warm IC (January Niño3 SST anomaly > 0.5°C), cold IC (January Niño3 SST anomaly < −0.5°C), and neutral IC (−0.5°C < January Niño3 SST anomaly < 0.5°C). Using step 2, we group the 100 January oceanic initial states into 24 warm cases, 29 cold cases, and 47 neutral cases. Here, it must be mentioned that the 0.5°C threshold corresponds to approximately half the SD (1.08°C) of the January Niño3 SST anomalies of the PI-CTL simulation of the IITM-ESM. The logic behind using the January SST anomalies in forming the three groups comes from the fact that the volcanic forcing ensemble runs are all started from January 1st ICs corresponding to the 100 identified cases. Further, we have verified that the January 1st Niño3 SST anomalies corresponding to the 24 warm, 29 cold, and 47 neutral cases are (>0.5°C), (<−0.5°C), and (−0.5°C < January Niño3 SST anomaly < 0.5°C), respectively (not shown).

The VRF1x experiment is a 100-member ensemble (i.e., 24 warm, 29 cold, and 49 neutral) simulation from 1 January 1883 through 31 December 1885 in which the model is forced by the Krakatoa VRF. Note that the eruption started in August 1883 and peaked in 1884 (Fig. 4A). The stratospheric AOD from CMIP6 (Supplementary Materials and Methods, section S2) for the Krakatoa eruption is used as forcing in the IITM-ESM.

We find that 41 of the 100 realizations in the VRF1x experiment transitioned to warmer-than-normal states (ΔNiño3 > 0.5°C), 25 members transitioned to cooler-than-normal states (Δ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 (fig. S7). It comes out that of the 41 realizations that transitioned to warmer-than-normal states in VRF1x, 8 cases started from warm ICs, 13 cases started from cold ICs, and 20 cases started from neutral ICs (fig. S7). By tracing back to the above 41 ICs, we repeated the 41-member ensemble experiment (VRF0.25x) with one-fourth VRF. Likewise, of the 25 realizations that transitioned to cooler-than-normal states (ΔNiño3 < −0.5°C) in VRF1x, we find that 8 cases started from warm ICs, 7 cases started from cold ICs, and 10 cases started from neutral ICs (fig. S7). We further repeated the 25-member ensemble simulations with quadrupled VRF (VRF4x).

Composite map of land surface temperature anomaly

Gridded paleoclimatic surface air temperature anomaly reconstructions for the last millennium (60) are used to construct the composite map of NH surface temperature anomalies for LVEs (i.e., magnitude of VRF > 3.7 W m−2) (12, 17). The composite anomaly is computed as follows:

Postvolcano = Average of NH gridded surface temperature anomalies for years 0 and +1 following a volcanic eruption.

Prevolcano = Average of NH gridded surface temperature anomalies for years −5, −4, −3, −2, and −1 before volcanic eruption.

The composite difference (postvolcano − prevolcano) map is shown in Fig. 1C.

SUPPLEMENTARY MATERIALS

Supplementary material for this article is available at http://advances.sciencemag.org/cgi/content/full/6/38/eaba8164/DC1

https://creativecommons.org/licenses/by-nc/4.0/

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is not for commercial advantage and provided the original work is properly cited.

REFERENCES AND NOTES

Acknowledgments: We thank the original data providers and the National Centers for Environmental Information (NCEI) for providing open access to paleoclimatic reconstructions via www.ncdc.noaa.gov/data-access/paleoclimatology-data. M. G. Yadava, H. Borgaonkar, and S. Feng provided the last-millennium paleoclimatic reconstructions for the Indian monsoon. In addition, we acknowledge the World Climate Research Programme (WCRP) for supporting the PMIP3 outputs and CMIP6 forcing datasets. We thank the two anonymous reviewers and the editor for providing important suggestions for improving the manuscript. Funding: This research is supported by the Ministry of Earth Sciences (MoES), Government of India. We thank the Director of the IITM for providing the necessary support. B.G. was partially funded by the DFG project IUCLiD (project no. DFG MA4759/8). R.V.D. has received financial support from the German Federal Ministry for Education and Research via the Belmont Forum/JPI Climate project GOTHAM (BMBF project no. 01LP1611A). R.K. and R.V. received funding support from the MoES, Govt. of India, via the Belmont Forum/JPI Climate Project GOTHAM to participate in the GOTHAM International Summer School, Potsdam, Germany, September 2017. This synergistic collaboration by J.K. and R.V.D. is within the scope of the IRTG 1740/TRP 2015/50122-0, funded by the DFG/FAPESP. Author contributions: M.S., R.K., and B.G. conceptualized and led the work; M.S., B.G., R.V.D., N.M., R.K., and J.K. contributed to data analysis (e.g., phase coherence, quantile-quantile analysis of ENSO-IM phase speeds, event synchronization, Bayesian analysis, and statistical testing) including validation and interpretation of results; M.S., A.D.C., P.S., A.G.P., and N.S. designed, performed, and analyzed the IITM-ESM targeted climate sensitivity experiments for the Krakatoa volcanic eruption; R.K., M.S., B.G., R.V., C.V., and J.K. wrote the manuscript; and all authors reviewed and edited the manuscript. Competing interests: The authors declare that they have no competing interests. Data and materials availability: All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. The code and data used to perform phase coherence analysis, event synchronization analysis, Bayesian analysis, generation of quantile-quantile plots of ENSO-IM phase speeds, and analysis of IITM-ESM model outputs are available at https://github.com/manmeet3591/fingerprint-volcano-enso-im. Any additional data or analysis codes related to this paper may be requested from the authors.
View Abstract

Stay Connected to Science Advances

Navigate This Article