We have mainly used a sub-set of the decadal hindcast outputs from the CMIP5 decadal runs (Meehl et al., 2009, 2014; Taylor et al., 2012; Smith et al., 2013; Boer et al., 2016). These runs have anthropogenic and natural forcing in them, and comprise several independent decade-long ensemble-runs by each model. For example, the four members of the first ensemble are initialized in 1960 and run up to 1970 (Figure 1). Such runs are available with initial conditions from 1960 to 2011. Together, we have 51 ensemble members of decadal hindcasts for each model, each differing from the other in the initialization year (from 1960 to 2011). For each model, we average multiple ensembles with same initialization dates. These averages are referred to as the model hindcasts for the initial conditions of that particular year. The models, whose hindcasts we consider, are the CanCM4 (Merryfield et al., 2013), MIROC5 (Watanabe et al., 2010), BCC-CSM (Wu et al., 2013), and the GFDL (Delworth et al., 2006; Zhang et al., 2007) models. The output data from all the models are divided into 10 groups depending on the year after initialization (Figure 1). Yr1, yr2, yr3… indicate 1, 2, 3 … year(s) after initialization.

We used Hadley Center Sea Ice and Sea Surface Temperature (HadISST; Rayner et al., 2003) and Ocean Reanalysis data from the European Center for Medium-Range Weather Forecasts (ORAS4; Balmaseda et al., 2013) for ocean temperatures. As for ENSO indices, we calculate the NINO3 and NINO4 indices as anomalies of SSTs in area averaged over the boxes 150–90°W, 5°S−5°N and 160°E−150°W, 5°S−5°N, respectively. We calculate the Indian Ocean Dipole Mode Index (IODMI) as the SST gradient between the western equatorial Indian Ocean (50–70°E and 10°S−10°N) and the south-eastern equatorial Indian Ocean (90–110°E and 10°S−0°). We take into consideration the IODMI only during the months of September to November as the IOD peaks during this season.

We have carried out correlation and regression analysis. We have used the two-tailed Student's t-test to determine the significance of the correlation coefficients. We ascertained the robustness of the skills by applying a boot-strapping test for 1,000 simulations (using NCL) and found our results to be significant at 90% for up to 4 years. For verifying the skill of persistence for both the models and observations, we use the auto-correlation (Supplementary Figure 1). To illustrate briefly, say in the case of CanCM4, we correlate the predicted IODMI at lead 1 with that at lead 2 to get the persistence skill at lead 1. Similarly, the persistence skill at lead 2 is obtained by correlating the predicted IODMI at lead 0 with that at lead 2, and so on. We have used the Empirical Orthogonal Function (EOF) analysis technique (von Storch and Zwiers, 1999) to determine the dominant patterns of spatio-temporal variability in the equatorial Indian Ocean, as an additional analysis to compare the simulated IODs with the observations. The EOF analysis is a multivariate statistical technique to calculate the eigenvalues and eigenvectors of a spatially weighted anomaly covariance matrix. The corresponding eigenvalues quantify the variance percent explained by each pattern, which is, by definition, orthogonal to the other patterns.

To have an estimate of the inter-annual variability of IOD that could be attributed to the Southern Ocean (20°W:50°E; 60°S:40°S), we projected the IODMI on to the corresponding area-averaged subsurface temperatures (vertically averaged over 300–800 m depth), i.e., through a regression analysis, and measured the goodness-of-fit (figure not shown).

ENSO events peak during the boreal winter (through December), while the IOD events peak during boreal fall (September–November). Our analysis shows that the hindcasts analyzed in this study (Figure 2A) predict the peak NINO3 index (Trenberth, 1997) significantly at leads up to 12–13 months, in agreement with earlier studies (Sun et al., 2018; Pal et al., 2020).

Surprisingly, the hindcasts from the MIROC5 and CanCM4 predict the IODMI during the fall season with significant lead prediction skill for at least 2–3 years (Figure 2B). The skill levels for most of the lead times are not only statistically significant at the 95% confidence level from a two-tailed Student's t-test but also are better than persistence (Supplementary Figure 1). While the lead prediction skills of the IODMI from the CanCM4 fall below 95% confidence level at 4–5-year leads, these skill levels are still significant at a 90% confidence level. The MIROC5 also shows a weak skill for the eighth year, still significant at 80% confidence level.

It is noteworthy that this skill is present despite the maximum lead predictability of only 1 year for the NINO3 index in all the models. We also ascertain that the maximum lead time skill for the NINO4 is similar to those for the NINO3 index. The ISM rainfall is an ocean–atmosphere coupled phenomenon occurring every year from June to September bringing copious amount of rainfall over the Indian subcontinent. ISM is strongly influenced by ENSO and IOD, along with other factors. Our study shows that ISM has no predictability in the analyzed models, presumably due to model inadequacies in representing monsoon processes (figure not shown). Nevertheless, rainfall predictions could be improved by exploiting known statistical relationships between ENSO-ISM or IOD-ISM (Jourdain et al., 2013; Swapna et al., 2015; Dutta and Maity, 2018).

To further ascertain whether the IODMI computed from the two hindcast datasets represents some of the observed features of the IOD (Yamagata et al., 2003), we compare the second leading simulated EOF modes of the SST from the hindcasts of MIROC5 and CanCM4 with that from the corresponding observations. Figure 3, derived from the EOF analysis at 1-year lead, shows that the models indeed capture the observed dominant IOD variance pattern associated with the EOF2 in terms of the location of the two centers of action. The simulated variance explained by this statistical mode from the MIROC5 is 12.2%, comparable with the corresponding observed value of 12% (Saji et al., 1999; Ashok et al., 2004). EOF2 of CanCM4 is 27.4%, which is higher than the expected value as the dipole-like structure shows more spread in the eastern Indian Ocean (Figure 3). The simulated EOF1, associated with Indian Ocean Basin mode, and the corresponding variance explained are also realistic.

A question arises as to what processes are responsible for the IOD predictability. Earlier studies have shown that the key component of improved decadal prediction lies in ocean physical processes, which internally generate decadal climate variability (Meehl et al., 2014; Farneti, 2017). Thus, the memory of the upper ocean (~800 m) provides an improved predictability of SST variability in models on decadal timescales (Yeager et al., 2012; Wang et al., 2014). Motivated by these, we looked for a source of predictability of IOD at several levels in the subsurface temperature (every 50 m). We found the maximum signal in the Southern Ocean at 300 to 800 m depth 7–10 years before the occurrence of the IOD event. From an analysis of the ORAS4 and SODA3 reanalysis datasets, we find a signal for the IODMI in the Southern Ocean. This is evidenced by significant correlations between depth-averaged 300–800 m temperatures in the Southern Ocean, which lead the IODMI by 10–6 years, respectively (Figure 4). The significant positive correlations in the Southern Ocean indicate that a positive temperature anomaly in the Southern Ocean leads to a positive IOD event after 8–10 years, whereas a similar lead negative correlation coefficient in the Southern Ocean indicates a negative IOD event after such a time period. The signal leading the IODMI first appears in the sub-surface temperatures of the Southern Ocean just south-west of Africa between the depths of 300–800 m, 10 years before the occurrence of the IOD event. This signal is seen to propagate toward the east along the Antarctic Circumpolar Current (ACC, or West wind drift) to about 50–60°E about 6 years prior to the IOD event. From here, the signal in the heat content further propagates north, along the east of Madagascar, all the way to the Somali coast in a pathway that is coincident with the East Madagascar undercurrent (EMUC)—an intermediate undercurrent at the depth of thermocline and below. The period of the propagation is in agreement with the EMUC as described by Nauw et al. (2008). The EMUC transports about 2.8 ± 1.4 Sv of intermediate water equatorward. The signal upwells to the surface near Somalia (Somali upwelling) along with its propagation north (Figure 5).

The propagation of the signal from the deeper layers in the extra-tropics to the surface layers in the equatorial region can be explained by the Indian Ocean Meridional Overturning Circulation (IOMOC; Wang et al., 2014). The time period of the propagation is in agreement with the IOMOC. From the Somali coast, the lead heat content signal is seen to move east, reaching the central equatorial Indian Ocean 3 years before the IOD event. The signal remains in the equatorial Indian Ocean for the last 2–3 years before the occurrence of the IOD event (Figure 6). Furthermore, our correlation in Figure 6 suggests an anomalous SST structure of opposite polarity in the equatorial Indian Ocean about a year before the occurrence of an IOD event. This is in conformation with Saji et al. (1999), who not only note a tight coupling between the intensity of the SST dipole anomaly and zonal wind anomaly but emphasize on a biennial tendency. The biennial tendency of the IOD is also noted by Meehl and Arblaster (2001), Ashok et al. (2003), etc. It is intriguing why the signal from the sub-surface at the equator does not manifest as an IOD event immediately. There are several studies that suggest potential drivers that can force an IOD such as the Indian summer monsoon, ENSO, and can account for the lead signals associated with the IOD with 1–2-year lead (Ashok et al., 2001; Ashok and Saji, 2007). As the current models do not show any lead skills of the Indian summer monsoon beyond a season, we can rule out that the long lead skills for the IODMI in the models come from the lead skills of monsoon or those related to ENSO (see Figure 2A). This might be explained by the fact that some IOD events are triggered by an atmospheric signal substantial enough to trigger the coupled evolution of the IOD (Shinoda et al., 2004; Saji, 2018) and hence this delay of 2–3 years.

A linear regression analysis carried out using HadISST and ORAS4 ocean temperature datasets (slope value = 0.46) indicates that the Southern Ocean signals at decadal lead explain about 18% of the inter-annual variability of IOD as suggested by a goodness-of-fit measure. To put it in perspective, ENSO, known as the most prominent driver of the Indian summer monsoon interannual variability, explains about 30% of the latter.

The source of decadal predictability for IOD in the MIROC5 and CanCM4 also apparently comes from the Southern Ocean (Supplementary Figure 3). We find statistically significant anomaly correlations between the depth averaged from surface until ~800 m hindcast heat content with the IODMI with the heat content leading at 10–5 years, conforming to our results from the reanalysis. The source of signal is the Southern Ocean for both the models, but the simulated signal path from here to the equatorial Indian Ocean are, however, weaker and remain unclear specifically in the MIROC5 model. However, the signals re-emerge finally 2–3 years before the occurrence of the simulated IOD. The models also do not capture the biennial tendency of the IOD as in observations. Further understanding of other possible factors affecting IOD variability on decadal scales is needed.