Global-scale interdecadal variability a skillful predictor at decadal-to-multidecadal timescales for Sahelian and Indian Monsoon Rainfall

In the present study, a sea surface temperature-based index named global-scale interdecadal variability (GIV) encompassing the combined variability of Atlantic multidecadal oscillation (AMO) and interdecadal Pacific oscillation (IPO) has been proposed. The warm phase of GIV exhibits a “cold AMO-like” pattern in the Atlantic basin and a “warm IPO-like” pattern in the Pacific basin. About 84% (R ~−0.914) of Sahelian and 42% (R ~−0.647) of Indian rainfall’s temporal variance is attributed to GIV, showing substantial improvement compared to the variance explained by AMO and IPO individually. The physical mechanism for GIV-rainfall teleconnection is related to a modification of the Walker circulation. Although there is a substantial degree of uncertainty in the current generation of state-of-the-art climate models from the Coupled Model Intercomparison Project Phase 5 (CMIP5), some still replicate the observed GIV’s spatial structure, its teleconnection, and associated physical mechanism. The results presented herein advance our knowledge about rainfall’s interdecadal variability and have imperative ramifications for developing skillful decadal predictions.


INTRODUCTION
Besides having variations on intraseasonal-to-interannual timescales, the Indian monsoon [1][2][3][4][5][6] and Sahelian monsoon 7-11 also exhibit pronounced variability at interdecadal (i.e., decadal-tomultidecadal; D2M) timescales, having alternate epochs of aboveand below-normal rainfall. Since the long-term floods and droughts are associated with the positive and negative interdecadal phases of rainfall, therefore, skillful prediction and a deeper understanding of these interdecadal swings could be of great value for policymakers and agricultural planning sector. However, prediction of these epochal changes in rainfall would require a better understanding of the causes of variability.
The observed interdecadal variability of rainfall over India and Sahel is considered as a result of large-scale forcing induced by the interplay of different oceanic basins 4,[12][13][14][15][16][17] . Hence, oceans play a vital role in the decadal variability and predictability of our climate because of their thermal and dynamical inertia. Meehl et al. 18 stated that several phenomena in the Pacific and the Atlantic Oceans could provide predictive skill on decadal timescales if initialized properly in models. In the Atlantic Ocean, the phenomenon is referred to as the Atlantic multidecadal oscillation (AMO), and in the Pacific, it is termed as the interdecadal Pacific oscillation (IPO).
The AMO, named by Kerr 19 and was first recognized in the mid-1990s, is a basin-scale variability of sea surface temperature (SST) in the North Atlantic, having a multidecadal period of 65-70 years, and an amplitude of 0.4°C [20][21][22] . The literature provides sufficient evidence that the warm phase of AMO is linked with critical climatic signals, for instance, the increased hurricane activity over the Atlantic 16,23,24 , the North American as well as European summer climate 25 , intensification of summer monsoon rainfall over Sahel 16,17,24,26,27 and India 4,12,13,15,16,28 . Likewise, the IPO is a basin-wide pattern of the main mode of SST variability in the Pacific at D2M timescales [29][30][31][32] , i.e., the Pacific-wide manifestation of the Pacific Decadal Oscillation (PDO) 33 , a North Pacific component of IPO. Many preceding studies demonstrated that the IPO also influences regional climate, for example, the warm phase of IPO is associated with wet periods over West and Central United States 34 and drought-like conditions over the Sahel 17,35 as well as over India 4,5,13,14 .
Recent studies highlight the integrated effect of these oceanic indices (i.e., AMO and IPO) on precipitation variability over conterminous United States 36 , Sahel 17 , and India 4 . Mohino et al. 17 reported that IPO's warm phase leads to drought-like conditions over the Sahel, while the AMO's warm phase enhances Sahelian rainfall. Joshi and Rai 4 revealed that the cold phase of IPO and the warm phase of AMO causes enhancement of rainfall over the Indian region. The studies discussed above insinuate a possible role of these oceanic indices on D2M rainfall variability over the Sahel and India. For the regions whose precipitation is governed by these indices, it is required to see their combined influence for better D2M rainfall prediction. So inspired by the previous studies 4,17,36 , herein, an index is proposed that encompasses their combined variability and will serve as a skillful predictor compared to AMO and IPO. Hence, this index will be of great importance and considered a worthwhile pursuit to understand the D2M variability and prediction of rainfall over respective regions.
In the present study, an attempt has been made to offer an SST-based index that contains D2M variability from the global oceanic basins, precisely the Pacific and North Atlantic Oceans. Further, it has been explored how this index is associated with the interdecadal variability of rainfall over Sahel and India, and what is the probable candidate mechanism responsible for such teleconnections.

Global-scale interdecadal variability
Using observational data, the present study suggests that a global SST illustrates a dominant mode of variability at interdecadal, i.e., D2M timescales (hereafter referred to as global-scale interdecadal variability (GIV)). The spatial structure and time series of GIV is obtained by doing empirical orthogonal function (EOF) analysis of the smoothed (11-year centered running mean) annual mean SST anomalies (SSTAs) over the region (40°S-60°N, 0°E-360°E) for the period 1906-2013. Prior to doing EOF analysis, the SST data are not de-trended; therefore, the first EOF (EOF-1) represents the global warming mode (Supplementary Fig. 1a) and its corresponding temporal coefficients (PC-1) divulges the global warming trends (i.e., ubiquitous warming over the oceans; Supplementary  Fig. 1b). The second EOF (EOF-2; Fig. 1a) illustrates the GIV mode and its corresponding PC (PC-2; Fig. 1b) shows how GIV evolves with time.
The warm phase of GIV (Fig. 1a) is typified as a "cold AMO-like" pattern in the Atlantic basin (characterized as a dipolar pattern, having colder SSTAs in the North Atlantic and warmer in the south; pattern just opposite to an AMO's warm phase, see Supplementary  Fig. 2a) and "warm IPO-like" pattern in the Pacific basin (i.e., having warmer SSTAs in the tropical Pacific, anomalies of opposite polarity in its western part, poleward of 25°and more eminent in the Northern Hemisphere; see Supplementary Fig. 3a).
The spatial correlation between the warm phase of GIV in the Atlantic (0°− 60°N, 75°W-5°W) and Pacific (40°S-60°N, 140°E-80°W ) basins with the warm phase of AMO and IPO is −0.935 and 0.831, respectively. The GIV index (Fig. 1b) reveals two warm (1901-1919 and 1971-1999) and two cold (1920-1970 and from 2000 onwards) phases during the analysis period. It is noteworthy that the exact locale of phase swapping may vary marginally, depending on the smoothing techniques and time windows one uses.
Likewise, the GIV index ( Fig. 1b), also illustrates an out-of-phase and in-phase relationship with AMO ( Supplementary Fig. 2b) and IPO ( Supplementary Fig. 3b) having a significant correlation −0.756 and 0.682, respectively. This implies that the fraction or proportion of variance of GIV explained by the AMO and IPO on an individual basis is 57.2% and 46.5%, respectively.
GIV epitomizes the combined interplay of AMO and IPO Now the question arises: can the GIV index be used to represent the combined variability of AMO and IPO? To answer this question, the multiple linear regression (MLR) technique is applied to see what percent of the variance of the GIV index is explained by the AMO and IPO?
But before applying the MLR technique, it is essential to verify that the two predictors (i.e., AMO and IPO) are independent and there is no effect of collinearity between them. The correlation between the two is very weak and insignificant (−0.179), suggesting that the AMO and IPO are independent and not strongly related to each other.
Further, to check the collinearity effect (for details see Methods), the AMO is regressed onto the IPO. The obtained values of coefficient of determination and variance inflation factor (VIF) are 0.0321 and 1.0332 (<10; critical threshold of VIF), respectively, signifying that there is no multicollinearity between the two independent variables. In general, multicollinearity requires strong inter-correlation between predictors and not just the non-zero inter-correlation. Now the GIV index is regressed concurrently onto the AMO and IPO for the analysis period (1906-2013) and the regression coefficients along with the adjusted coefficient of determination (R 2 adj ), statistics compensating the number of variables in the model and upsurges only if the auxiliary variables contribute significantly to the model, are computed.
The combined effect of AMO and IPO as shown by the MLR model is −0.364*AMO + 0.566*IPO that accounts for 87.7% of the variance of GIV (having a multiple correlation coefficient, R = 0.937) suggesting that the GIV index can be used as a representative of the combined variability of AMO and IPO. The negative and positive sign of AMO and IPO's regression coefficients further confirm their out-of-phase and in-phase relationship with GIV, as discussed above.

GIV-rainfall teleconnection
Now to examine the role of GIV in modulating the rainfall, the June-September (JJAS) seasonal precipitation anomalies are regressed onto the standardized GIV index (see Fig. 2a). The GIV-precipitation regression pattern reveals that during the cold phase of GIV (see Methods), the Sahel and central India experiences above-normal rainfall. On close observation, one can visually see that the GIV-precipitation regression pattern ( Fig. 2a) shows slightly higher precipitation anomalies over the respective regions as compared to the individual AMO-and IPOprecipitation regression patterns (Fig. 2b, c). Now to explore this quantitatively, the temporal variability of smoothed (11-year centered running mean) precipitation areaaveraged over the Sahel (7.5°N-20°N, 17°W-40°E; area enclosed within the red dashed-line in Fig. 2a) and central India (15°N-27.5°N , 70°E-87.5°E; area enclosed within the red solid-line in Fig. 2a) are plotted along with GIV, AMO, and IPO indices (see Fig. 2d, e) to scrutinize the proportion of variance in the rainfall over the respective regions that is explained by the GIV, AMO, and IPO on an individual basis as well as by the AMO and IPO on a combined basis.
It is perceived that during the cold (warm) phase of GIV, the Sahel and central India have received above-(below-) normal rainfall, which signifies the out-of-phase relationship between GIV and rainfall over respective regions, since high significant negative correlations −0.914 and −0.647 between the two are observed. As compared to the AMO and IPO, the rainfall over the respective regions (especially Sahel) closely follows the GIV, which is well depicted from the strength of the correlations (see Fig. 2d, e).
Noteworthy, none of the two predictors (i.e., AMO and IPO) by themselves can account for more than 42.4% and 25.8% of the rainfall variance for the Sahel and central India, respectively; whereas on considering both the AMO and IPO as predictors, the proportion of the variance of rainfall accounted for by regression increases to 69% and 41.7% for the respective regions. However, if one uses the GIV index as a predictor, then the percent of the variance of rainfall explained by it shows substantial improvement in the case of Sahel (i.e., 83.6%) and quite similar results in the case of central India (i.e., 41.9%).
Consequently, the regions (especially the ones considered herein) whose precipitation are governed by the AMO and IPO both, then for those regions, it is essential to see the combined influence of these oceanic indices for better D2M prediction of rainfall, and as discussed above, the GIV index is a better representation of these oceanic indices together.

Large-scale circulation features and physical mechanism
In principle, the individual contributions of these oceanic indices to the African and the Indian monsoons have been discussed in detail in many previous studies. To identify the combined effect of these oceanic indices and to scrutinize the atmospheric circulation patterns, as well as the physical mechanism allied with GIV, seasonal (JJAS) mean anomalies of atmospheric variables from NOAA-CIRES 20CR are regressed onto the standardized GIV index for the period 1906-2013.
The observed GIV winds regression pattern at a lower level (vector; Fig. 3a) divulges strong inflow in West Africa from the Atlantic and the Guinea coast in the form of Westerly African Jet (WAJ) and in the Indian sub-continent from the Arabian Sea in the form of Somali Jet (i.e., the cross-equatorial flow). These jets strongly influence the West African and Indian monsoons by bringing abundant moisture from the warm waters of the Atlantic and the tropical Indian Ocean, respectively. This monsoonal inflow is mainly induced by the relative local land-sea temperature contrast.
Both the African and the Indian monsoon systems are typified by seasonal atmospheric circulation, which is rather complicated in the West African region where two easterly jets, one at the midtropospheric level and the other at the upper level, are involved 37 .
At the mid-tropospheric level, the African easterly jet (AEJ), driven by an anticyclonic circulation right above the Saharan heat  low, is amongst one of the crucial parameters affecting the West African monsoon (see Supplementary Fig. 4). As reported in many previous studies [38][39][40][41][42][43][44][45] , the latitudinal position of AEJ (mainly governed by the low-level WAJ) rather than its intensity plays a pivotal role in controlling the precipitation over the Sahel. The GIV winds regression pattern at 600 hPa (vector; Fig. 3b) illustrates weaker than the normal mid-tropospheric easterlies south and under the main axis of the AEJ. On close observation, one can see that the AEJ associated with anticyclonic circulation seems to be displaced northeastward, which enhances the vertical wind shear over the Sahel that in turn strengthens the baroclinic instability, probably the dominant mechanism over the respective region.
The Tropical Easterly Jet (TEJ), discovered by Koteswaram 37 , is an inter-regional circulation feature linking the Sahelian and Indian monsoons and is observed unfailingly in the upper tropospheric levels (100-200 hPa) during the JJAS season. The GIV winds regression pattern at 150 hPa (vector; Fig. 3c) reveals strong easterlies extending from Southeast Asia across the Indian Ocean and Africa to the Atlantic, i.e., nearly over half of the Earth's circumference.
To see the partial impacts of AMO and IPO on the large-scale circulation features, the horizontal winds at respective levels are regressed onto the standardized AMO and IPO indices, keeping IPO and AMO constant (see Supplementary Figs. 5 and 6). The large-scale circulation features allied with AMO and IPO also reveal the WAJ and cross-equatorial flow at the lower level and easterlies at the upper level but weaker in magnitude than the GIV related patterns (see Fig. 3).
Comparing GIV wind regression patterns at lower and upper levels, one can notice a baroclinic type response over the Sahel and Indian Ocean (notably in the Somali Jet region). To confirm this, the GIV zonal wind shear (u850-u150) regression is computed (see Supplementary Fig. 7) that depicts easterly vertical wind shear over India and in the vicinity of the Western Indian Ocean extending up to the Atlantic through West Africa, which is associated with above-normal rainfall over India and Sahel, respectively (see GIV-rainfall regression pattern Fig. 2a).
Further, to better understand the circulation features connecting Sahel and India with the proximate regions, the velocity potential, a measure of divergent flow and is often used as a proxy for large-scale ascent and descent, is calculated. At the upper level, a region of positive potential has converging winds, which exemplify strong subsidence, i.e., sinking motion at a lower level and vice versa.
The GIV velocity potential regression pattern illustrates that the cold phase of GIV is allied with anomalous convergence over the central-eastern tropical Pacific and divergence over the warm-pool region as well as over the extended Indian monsoon region and the Sahel at upper levels (Fig. 4b) and with anomalous divergence and convergence over the respective regions at lower levels ( Fig. 4a). At a lower level, the anomalous divergence (i.e., subsidence) over the central-eastern tropical Pacific connotes anomalous Walker cell that in turn fortifies the zonal overturning circulation. Because of this subsidence, the atmosphere is highly stable over the central-eastern tropical Pacific, which is unfavorable to and limits the occurrence of deep convective clouds and precipitation to occur. In contrast to the central-eastern tropical Pacific, the atmosphere is highly unstable over the warm pool and its neighboring regions, which favors the formation of deep convective clouds and rainfall to occur frequently. This anomalous circulation reinforces the easterly (westerly) and westerly (easterly) trade winds across the tropical Pacific (Atlantic) in the lower and upper atmosphere, respectively.
Further, we have also plotted the partial regression of the JJAS anomaly of velocity potential at respective levels onto the standardized AMO and IPO indices, keeping IPO and AMO constant (see Supplementary Figs. 8 and 9). The AMO and IPOrelated velocity potential regressions also reveal similar patterns but of slightly weaker magnitude, and the centers of positive and negative potentials in the case of AMO seem to be marginally displaced compared to GIV (see Fig. 4).

Model's fidelity in simulating GIV and GIV-rainfall teleconnection
To examine the fidelity of the state-of-the-art climate models in replicating the observed GIV's spatial structure, its teleconnection with rainfall over the regions considered, and the allied mechanism, the historical simulations of 30 models (i.e., the atmosphere-ocean global climate models (AOGCMs) and the Earth system models (ESMs)) from the Coupled Model Intercomparison Project Phase 5 (CMIP5) are used for the period 1901-2004. The GIV index for forced simulations is derived using the same methodology as discussed earlier. Since an 11-year centered running mean is applied to separate the interannual signal from interdecadal; therefore, the analysis discussed herein is limited for 1906-1999. Figure 5 elucidates that out of 30 CMIP5 models under consideration, very few can reproduce the observed GIV's spatial structure with a correlation greater than 0.25 and a normalized standard deviation between 0.95 and 1.05. The following models, namely, CanESM2, HadGEM2-AO, IPSL-CM5A-MR, MIROC-ESM, MIROC-ESM-CHEM, and MRI-ESM1, fulfill the above criterion. Hence, the multi-model ensemble (MME) mean of these six good models (MME good) is only considered for further analysis. However, the MME mean obtained by averaging across all models also divulges an explicit spatial structure of GIV having a pattern correlation of 0.35 with observation but smaller in magnitude because of averaging across independent model realizations.
To get further insight into the crucial elements in the modeled GIV pattern, the ensemble mean of second EOF across all good models (MME good; Fig. 6a) has been computed. Consistent with the observation, the ensemble mean of good models also exemplifies a cold phase of GIV (i.e., IPO negative phase and AMO positive phase). On close observation, it can be seen that even the ensemble mean of good models shows limitations in reproducing the precise comma-shaped structure of AMO as reported in earlier studies 15,46 , i.e., though it illustrates the warm SSTAs (slightly weaker in magnitude as compared to observations) over the sub-polar region, but very weak or non-existent warm SST loadings over the sub-tropical North Atlantic. Consistent with observation (see Fig. 2a), the ensemble mean of GIV-rainfall regression patterns across all good models (MME good; Fig. 6b) also divulges above-normal rainfall over Sahel and India during the cold phase of GIV, indicating the out-of-phase relationship between the two over the respective regions.
Further, to examine the fidelity of selected CMIP5 models in reproducing the GIV related atmospheric circulation patterns and physical mechanism, the seasonal (JJAS) mean wind anomalies are regressed onto the standardized GIV index. The ensemble mean of GIV winds regression patterns across all good models at the lower level ( Supplementary Fig. 10a) show WAJ entering West Africa from the Atlantic and the Guinea coast and the cross-equatorial flow in the Indian Ocean that seems to be shifted slightly northward and weaker in magnitude, compared to observations (Fig. 3a).
At the upper level, the ensemble mean of GIV winds regression patterns ( Supplementary Fig. 10b) also divulges strong easterlies extending from Southeast Asia to the Atlantic; however, the   (40°S-60°N, 0°E-360°E) across all good (MME good) CMIP5 models. Ensemble means of regression patterns, obtained by averaging the regression maps of JJAS seasonal anomaly of b precipitation (mm day −1 ), c velocity potential at 850 hPa (10 6 m 2 s −1 ), and d velocity potential at 150 hPa (10 6 m 2 s −1 ) onto the standardized GIV index across all good (MME good) CMIP5 models. The gray contours in Fig. 6b indicate the regions where the regression coefficient is statistically significant at the 90% confidence level. The vectors in Fig. 6c, d represent the divergent wind (m s −1 ). The units of regression coefficients are in per standard deviation.
pattern over the Indian sub-continent seems to be shifted slightly northwards and weaker in magnitude as compared to observations (Fig. 3c).
Consistent with observation, the ensemble mean of GIV velocity potential regression patterns for good models exemplify positive potential over the central-eastern tropical Pacific and negative over the warm-pool region as well as over the vast Indian monsoon region at upper levels (Fig. 6d) and a reverse pattern at lower levels (Fig. 6c).
Thus, the CMIP5 models that reasonably reproduce GIV's spatial structure also simulate its teleconnection with rainfall over the regions considered and allied mechanism.

DISCUSSION
In the present study, an SST-based index, referred to as GIV, representing the combined variability of AMO and IPO, has been proposed. Further, it has been explored how the proposed index is associated with the Sahelian and Indian monsoon rainfall and what is the probable mechanism responsible for such teleconnections.
The GIV's warm phase epitomizes a "cold AMO-like" pattern in the Atlantic basin and a "warm IPO-like" pattern in the Pacific basin. The results reveal that 87.7% of GIV's variance is explained by the AMO and IPO, suggesting that the GIV index can be used as a representative of their combined variability. The interdecadal variability of rainfall over the regions considered herein closely follows GIV in an opposite manner. This study bestows GIV as a better predictor, explaining 83.6% and 41.9% of the Sahelian and Indian rainfall variance that illustrates considerable improvement or added value compared to the variance explained by the AMO and IPO individually or on a combined basis.
The large-scale circulation features reveal that during the cold phase of GIV, the WAJ becomes more intense that drives more moisture flux inside the African continent, the TEJ becomes stronger, and the AEJ becomes weaker, enhancing the vertical wind shear over the Sahel that in turn strengthens the baroclinic instability and resulting in above-normal rainfall over it. Likewise, India will also have above-normal rainfall during GIV's cold phase, which is evident from large-scale circulation patterns associated with it, i.e., stronger cross-equatorial flow resulting in intense low-level westerlies that brings abundant moisture from the Indian Ocean and strong easterlies (i.e., TEJ) at the upper level. The velocity potential regression pattern reveals strong subsidence (convection) over the central-eastern (western) Pacific, which denotes anomalous Walker cell that strengthens the zonal overturning circulation. Due to this subsidence, the atmosphere is highly stable over the centraleastern tropical Pacific, limiting the formation of deep convective clouds and precipitation to occur. On the contrary, the atmosphere is highly unstable over the western Pacific and its neighboring regions, favoring the formation of deep convective clouds and rainfall to occur. Thus, the physical mechanism for GIV-rainfall teleconnection is related to a modification of the Walker circulation.
Moreover, we have also scrutinized the fidelity of 30 CMIP5 models (i.e., AOGCMs and ESMs) in simulating the observed GIV's spatial structure, its teleconnection with rainfall over the regions considered, and the allied mechanism. Out of 30 CMIP5 models under consideration, only six reasonably reproduce the GIV's spatial structure and simulate its teleconnection with rainfall over the Sahel and India. Nevertheless, these CMIP5 models also replicate the observed atmospheric circulation features and the physical mechanism allied with GIV. Thus, to better understand rainfall's interdecadal variability and decadal predictions over the regions considered, the models must simulate the GIV skilfully.
Decadal climate prediction, an emerging branch of climate science, aims to bridge the gap between short-term weather forecasts and long-term climate projections 47 and is widely recognized as a key for the applications to water resources, agriculture, energy, and infrastructure development 48 . For instance, the decadal variability influences the spatiotemporal patterns of rainfall that cause a problem for long-term agricultural planning; hence, in the agricultural sector, selecting future cultivation areas requires future climate information, which decadal predictions can provide. Likewise, the energy infrastructure (e.g., solar parks, wind farms, and dams) involves long-term investments in the energy sector, so decadal predictions can support decision-making on such investments to minimize losses. The infrastructure sector is also heavily affected by the climate variability and change that alters the frequency and intensity of extreme events and causes long-term changes in climatic conditions, affecting critical infrastructures like railways and roads as well as the power and water supply networks. Notably, changes in precipitation and hydrological conditions result in more extended drought periods that influence water availability in river catchments and dams. Improved decadal predictions significantly influence improved planning that suggests development strategies, resource allocation, and the set up of national priorities. So if the future depends on the weather, then decadal predictions will allow better planning.
This study aspires to improve our understanding of the processes driving the D2M climate variability and quantify how the interactions of these processes provide and/or limit predictability. Furthermore, the results conferred herein are vital for developing reliable decadal predictions at regional scales, which depends upon the better understanding of the linked mechanisms and identification of climate patterns that can potentially provide predictive skill several years in advance if initialized properly in models 18 . This will be beneficial for society to reduce the adverse impact and economic loss by increasing preparedness.

METHODS Data
To define an AMO and IPO indices as well as an index, proposed herein (i.e., GIV), representing the combined variability of the two, the monthly SST data at 2°× 2°resolution is obtained from the National Oceanic and Atmospheric Administration (NOAA) Extended Reconstructed SST, version 5 (ERSST.v5) 49 . Further, to investigate the impact of these oceanic indices on the Sahelian and Indian monsoon rainfall, the global gridded precipitation data (resolution 0.5°× 0.5°) over the land surface from Climatic Research Unit time series version 4.03 (CRU TS v. 4.03) 50 is used. Moreover, to scrutinize the large-scale circulation features and the probable physical mechanism allied with the proposed index, the monthly means of zonal and meridional winds (resolution 1°× 1°) at 850, 600, and 150 hPa are obtained from the NOAA Cooperative Institute for Research in Environmental Sciences (CIRES) Twentieth Century Reanalysis (20CR) V3 51 .
To scrutinize the fidelity of state-of-the-art climate models in reproducing the GIV's spatial structure, its teleconnection with rainfall, and the associated mechanism, the monthly means of variables like SST, precipitation, and horizontal winds from the historical simulations of 30 models that participated in the CMIP5 database are used (for details regarding modeling groups and resolution see Supplementary Table 1). The historical simulations are forced with the observed atmospheric composition changes that comprise both natural (e.g., solar forcing, volcanic impacts, aerosols, and emanations of short-lived species and their precursors) and anthropogenic components (i.e., anthropogenic aerosols and greenhouse gases) along with the time-varying estimates of land cover 52 . Although many models provide several realizations, for fair comparison as well as for an equal weightage, the first ensemble member (i.e., r1i1p1) of each CMIP5 model has been used in this analysis as followed by earlier studies 14,15,53,54 . These simulations are freely available at the website http://esgf-index1.ceda.ac.uk, maintained by Earth System Grid Federation (ESGF). The resolution of each dataset differs; therefore, for ease of comparison, the model outputs and the observational data (except the precipitation data from CRU) are interpolated into a common latitude-longitude grid (2.5°× 2.5°) by bilinear interpolation.

Methodology
The analysis presented herein is primarily based on regression maps, obtained by linearly regressing the field of interest onto the normalized index [53][54][55] . Thus, the units of regression maps and the field used for regressions will be the same. It is worth mentioning that all the fields are linearly de-trended before applying regression analysis so that the trends in the data do not affect the results. For computing the statistical significance of regression coefficients in the case of observations and reanalysis data, a two-tailed Student's t test is used; whereas in the case of MME regressions, a Student's one-sample t test is applied.
In the present study, IPO is defined as the spatial pattern (Supplementary Fig. 3a) and time series (Supplementary Fig. 3b) of the first EOF (EOF-1) of smoothed (11-year centered running mean) de-trended annual mean SSTAs computed over the Pacific basin (40°S-60°N, 140°E-80°W) for the period 1906-2013. Before doing EOF analysis, the annual mean SSTAs for the period 1901-2018 are linearly de-trended to eradicate the global component of the anthropogenic forcing.
Herein, the AMO index (shown in Supplementary Fig. 2b) is obtained by applying 11-year centered running mean on linearly de-trended annual mean SSTAs area-averaged over the North Atlantic basin (0°-60°N, 75°W-5°W ). The definition for constructing the AMO index is similar to those used in previous studies 25, [56][57][58] . The spatial pattern of AMO (shown in Supplementary Fig. 2a) is then constructed by regressing the annual mean SSTAs onto the standardized AMO index for the period 1906-2013.
For computing these indices (i.e., IPO, AMO, and GIV), an 11-year centered running mean is applied on annual SSTAs (1901-2018) to separate the interannual (referred to as high frequency) signal from interdecadal (referred to as low frequency). Thus, restricting the observed analysis from 1906 to 2013.
Before applying the MLR technique, it is vital to confirm that the predictors are independent and there is no effect of collinearity between them. To check the collinearity effect, the VIF, a statistic for identifying multicollinearity in a matrix of predictor variables, is computed using the following equation where R 2 j is the coefficient of determination obtained by regressing x j on the remaining k À 1 explanatory variables. The value of R 2 j will be close to one and VIF j will be large, if x j is linearly dependent on the other explanatory variables. If VIF exceeds 10, then there is strong evidence of collinearity in the set of explanatory variables.
It is to be noted that prior to doing regressions, the observed GIV index is multiplied by −1 to see the effect of the cold phase of GIV, which consists of AMO's warm phase and IPO's cold phase. This is conferred by regressing the observed GIV index multiplied by −1 onto the annual SSTAs (figure not shown). Similarly, the observed IPO index is also multiplied by −1 to see its cold phase effect.

DATA AVAILABILITY
The observational datasets used in this study are the monthly SST data from National Oceanic and Atmospheric Administration (NOAA) Extended Reconstructed SST version 5 (ERSST.v5) 49 at 2°× 2°resolution and the global gridded precipitation data over the land surface from Climatic Research Unit time series version 4.03 (CRU TS v. 4.03) 50 at 0.5°× 0.5°resolution. The monthly means of zonal and meridional winds at 850, 600, and 150 hPa are obtained from the NOAA Cooperative Institute for Research in Environmental Sciences (CIRES) Twentieth Century Reanalysis (20CR) V3 51 at 1°× 1°r esolution. The historical CMIP5 simulations used herein are freely available at the website http://esgf-index1.ceda.ac.uk maintained by ESGF.