Emerging Skill in Multi-Year Prediction of the Indian Ocean Dipole

The Indian Ocean Dipole is a leading phenomenon of climate variability in the tropics, which affects the global climate. However, the best lead prediction skill for the Indian Ocean Dipole, until recently, has been limited to ~6 months before the occurrence of the event. Here, we show that multi-year prediction has made considerable advancement such that, for the first time, two general circulation models have significant prediction skills for the Indian Ocean Dipole for at least 2 years after initialization. This skill is present despite ENSO having a lead prediction skill of only 1 year. Our analysis of observed/reanalyzed ocean datasets shows that the source of this multi-year predictability lies in sub-surface signals that propagate from the Southern Ocean into the Indian Ocean. Prediction skill for a prominent climate driver like the Indian Ocean Dipole has wide-ranging benefits for climate science and society.


INTRODUCTION
The Indian Ocean Dipole (IOD) is an inter-annual coupled ocean-atmosphere phenomenon in the tropical Indian Ocean that peaks during the boreal fall season (September to November; SON). Positive IODs are associated with anomalously warmer western Indian Ocean and anomalously cooler eastern Indian Ocean. Negative IOD events are associated with opposite anomalous sea surface temperatures (SSTs) across the Indian Ocean (Saji et al., 1999;Webster et al., 1999;Murtugudde et al., 2000). The IOD is a well-known driver of global climate (Saji and Yamagata, 2003). In addition to affecting the neighboring Maritime Continent to the east and East Africa to the west (Behera et al., 2005), the positive IOD events, for example, have been associated with reduced rainfall over western and southern Australia Ashok and Saji, 2007;Cai et al., 2011), enhanced seasonal Indian summer monsoon (ISM) rainfall (Ashok et al., 2001;Ashok and Saji, 2007), and climates of even distant regions of South America (Chan et al., 2008) and Europe (Hardiman et al., 2020). In addition to its own coupled dynamics (Gualdi et al., 2003;Yamagata et al., 2004;Behera et al., 2006;Ha et al., 2017;Tanizaki et al., 2017;Saji, 2018;Marathe et al., 2021), the IOD is suggested to be triggered by other inter-annual processes such as the El Niño-Southern Oscillation (ENSO), the dominant climate driver from the tropics.
Attempts to predict the IOD on a seasonal scale have been on-going for about two decades (Iizuka et al., 2000;Shinoda et al., 2004;Luo et al., 2007;Doi et al., 2016) and have met with relatively shorter lead time skills, the highest being of 4 months. This is in contrast to the 12-17 months of lead prediction skill for the ENSO as shown in recent studies (Barnston et al., 2012;Park et al., 2018;Tang et al., 2018). Several papers claim that prediction skills for the IOD and ENSO are linked and intertwined (Wajsowicz, 2005;Izumo et al., 2010;Luo et al., 2010).
The emerging discipline of decadal prediction, i.e., prediction of climate information for the near future, shows great promise for societal needs and economic policymakers. Decadal prediction lies between weather prediction and climate projection. Therefore, decadal prediction has to resolve initial value problems like in weather prediction and seasonal to inter-annual forecasts and climate variability and trends effective from boundary conditions, including slow varying oceanic processes, snow cover, and anticipated changes in anthropogenic greenhouse gases and aerosols (Meehl et al., 2009(Meehl et al., , 2014Smith et al., 2013;Boer et al., 2016). Interestingly, the inter-annual ENSO and the IOD can also be modulated by decadal processes (Ashok et al., 2004;Tozuka et al., 2007). Such processes should enhance the decadal prediction skill of these inter-annual events.
Given the limited lead skill of the IOD on seasonal scales, its prediction at multi-year lead times is exceptionally challenging. Remarkably, in this study, we show that two model hindcasts from the CMIP5 (Coupled Model Intercomparison Project-Phase 5) decadal prediction datasets are found to show successful prediction skill for IOD events. These sets give us an opportunity to track the lead prediction skill of an event at a lead of up to 2 years.

Data
We have mainly used a sub-set of the decadal hindcast outputs from the CMIP5 decadal runs (Meehl et al., 2009(Meehl et al., , 2014Taylor 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.

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

Lead Prediction Skills for ENSO and IOD
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 FIGURE 1 | The output data from all the models are divided into 10 groups depending on the year after initialization. Yr1, yr2, yr3… indicates 1, 2, 3 … year(s) after initialization, which is represented by individual colors.
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).

Fidelity of the Simulated IOD
To further ascertain whether the IODMI computed from the two hindcast datasets represents some of the observed features of the IOD , 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.

Source of the IOD Prediction Skill
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 depthaveraged 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  Frontiers in Climate | www.frontiersin.org 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 FIGURE 4 | Spatial distribution of the anomaly correlation of HadISST derived IODMI with vertically averaged subsurface temperatures from ORAS4 (over 300-800 m depth) from previous years averaged annually. The "lag number" is the years over which the IOD lags the Southern Ocean signal. In all panels, significant correlations at 95% confidence are shown as contour lines (value at 0.27). 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.

DISCUSSION
Our analysis of the CMIP5 decadal retrospective prediction products for the 1960-2011 period shows that two CMIP5 decadal prediction systems exhibit multi-year skill in predicting FIGURE 5 | Vertical distribution of anomaly correlation of IODMI with previous years' subsurface temperatures, averaged annually over longitudes 50-100 • E. Significant correlations at 90% confidence (0.23) and at 95% confidence (0.27) are contoured. the IODMI, higher than the 17-month lead skills for the ENSO, a leading climate driver. The physical processes responsible for the skill levels seems to originate in the Southern Ocean up to a decade before emerging in the Indian Ocean, suggesting that this multi-year skill could be extended even further. To our knowledge, this is the first study to show decadal skill in predicting the IOD.
Skill levels for surface fields and, in particular, rainfall, are not statistically indistinguishable from zero. This is in common with long-term predictions of other climate phenomena (Scaife and Smith, 2018) and indicates the relative immaturity of this form of forecasting. It remains a challenge to both improve models and to fully exploit decadal prediction skill. Nevertheless, the potential for long-term planning in both public and private sector organizations that could be facilitated by such skillful predictions is great.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding authors.

AUTHOR CONTRIBUTIONS
FF, KA, and MC formulated the hypothesis and co-wrote the manuscript. FF conducted the analyses. All authors reviewed the final version of the manuscript.

FUNDING
FF thanks MC and NERC/MoES SAPRISE project (NE/I022841/1) for the travel grant to visit MetOffice, UK. FF acknowledges CSIR for the JRF & SRF fellowship. FF also thanks C T Tejavath for assistance in data processing.