Information Entropy as Quantifier of Potential Predictability in the Tropical Indo-Pacific Basin

Global warming is posed to modify the modes of variability that control much of the climate predictability at seasonal to interannual scales. The quantification of changes in climate predictability over any given amount of time, however, remains challenging. Here we build upon recent advances in non-linear dynamical systems theory and introduce the climate community to an information entropy quantifier based on recurrence. The entropy, or complexity of a system is associated with microstates that recur over time in the time-series that define the system, and therefore to its predictability potential. A computationally fast method to evaluate the entropy is applied to the investigation of the information entropy of sea surface temperature in the tropical Pacific and Indian Oceans, focusing on boreal fall. In this season the predictability of the basins is controlled by two regularly varying non-linear oscillations, the El Niño-Southern Oscillation and the Indian Ocean Dipole. We compute and compare the entropy in simulations from the CMIP5 catalog from the historical period and RCP8.5 scenario, and in reanalysis datasets. Discrepancies are found between the models and the reanalysis, and no robust changes in predictability can be identified in future projections. The Indian Ocean and the equatorial Pacific emerge as troublesome areas where the modeled entropy differs the most from that of the reanalysis in many models. A brief investigation of the source of the bias points to a poor representation of the ocean mean state and basins' connectivity at the Indonesian Throughflow.


INTRODUCTION
In the past decade, our theoretical understanding of the physics of the climate system has advanced in fundamental ways. These advancements proceeded in parallel with model improvements and computing capabilities. Understanding and especially predicting climate change at regional or local scales-the scales that are relevant to society-remains, however, challenging: regional climate change is influenced by the large-scale climate and, at the same time, feeds back to the global scale but these interactions are poorly represented in models.
In this work, we focus on the tropical Pacific and Indian Oceans, where the El Niño Southern Oscillation (ENSO; Bjerknes, 1966Bjerknes, , 1969 and the Indian Ocean Dipole (IOD; Saji et al., 1999;Webster et al., 1999) control the largest portions of the variance at interannual scale (Figure 1). They impact key variables of societal relevance, from surface temperature to precipitation and the frequency of extreme events such as tropical cyclones, typhoons and droughts. These climate modes are of the uttermost importance not only for water and food security, but also for global health, as they modulate, for example, malaria occurrences in India (Dhiman and Sarkar, 2017;Anyamba et al., 2019), Indonesia (Kovats, 2000) and Africa (Hashizume et al., 2012;Kreppel et al., 2019). Assessing their potential predictability and quantifying how this predictability is simulated in climate models and projected to change in the future, is a priority.
ENSO and IOD play a crucial role in global climate variability. In light of this importance, their interactions have been the subject of numerous studies. A recent review of ENSO teleconnections can be found in Yeh et al. (2018). The IOD develops in late summer with a November maximum. Since its discovery, observational and modeling studies have focused on its linkages with the Asian Summer Monsoon (Behera et al., 1999;Ashok et al., 2001;Guan et al., 2003;Saji and Yamagata, 2003), its modulation of the Indian Summer Monsoon and ENSO (Ashok et al., 2001); and its teleconnections over North America, Australia, South-Africa (Li and Mu, 2001) and East Africa (Black et al., 2003). The regional and global influences of these modes cannot be overemphasized, yet a realistic simulation of the characteristics and teleconnections of these modes remains challenging in state-of-the-art climate models (Weller and Cai, 2013), precluding the reliability of future projections.
Here, we introduce the climate science community to a quantifier of complexity, or entropy, recently proposed in the literature that allows for distinguishing regular, chaotic, and random behaviors in time-series. It provides a quick framework to investigate potential predictability and verify how well climate models represent it. We test it using an observationally-based reanalysis product and model data-sets, we verify its robustness and finally use it to exemplify the challenges implicit in using global climate models-here the integrations from the Coupled Model Intercomparison Project phase 5, CMIP5- (Taylor et al., 2012) to investigate present and future predictability. We focus on the Indo-Pacific Ocean in boreal fall, when the IOD variance is maximized and many countries are impacted by its and ENSO modulation of rainfall and temperatures over the surrounding land masses.
The rest of this paper is organized as follows: section Data provides information about the observational data and models used; section The Entropy Quantifier: Recurrence Plots and Information Entropy focuses the information entropy quantifier; sections Entropy in Reanalysis Data, CMIP5 Historical Runs and CMIP5 Projections discusses the entropy of historical and future projection focusing on the boreal fall season, and section Entropy and the ENSO-IOD Relation analyzes possible sources of model bias in the historical period. A summary of the findings concludes the work.

DATA
We analyze monthly mean outputs from 15 CMIP5 models and consider both historical runs and future projections following the representative concentration pathway (RCP) 8.5 scenario. For each model, we use at least 1 member for the historical analysis and its evolution in the RCP case. For 10 models we verified that internal model variability does not impact the outcome of our analysis by exploring 3 (or 2 when three runs were not readily available) members in the historical period. The output variable we focus on is sea surface temperature (SST), and further consider sea level pressure (SLP) and thermocline depth (Z20) in investigating model divergence.
In terms of observational datasets, we use the monthly mean subsurface temperature from the Simple Ocean Data Assimilation (SODA) version 3.4.2, which was forced by the European Center for Medium-range Weather Forecasts (ECMWF) re-analysis ERA-Interim (Carton et al., 2018). SODA has a spatial resolution of 0.25 • x 0.25 • with 50 vertical levels, and spans the 1980-2017 period. For cross validation, we also considered fields from the Ocean Re-Analysis System 4 (ORAS4) (Zuo et al., 2019). Results obtained using SODA and ORAS4 are nearly identical and consistent. For brevity, we focus our discussion on SODA. Sea level pressure is from the monthly mean ERA-Interim (Dee et al., 2011). All data have been detrended and re-gridded into a uniform 1 • longitude x 1 • latitude spatial resolution before analysis. For the observational datasets, we look at the period 1980-2018 for temporal consistency among them, while in the CMIP5 models we consider the last 39 years of the historical period -same period length as in In the following, the IOD strength and variability is quantified by the IOD Index, which is the difference in the SST monthly anomalies (SSTAs) averaged over the western tropical Indian Ocean (WTIO) (50 • -70 • E, 10 • S−10 • N) and in the SETIO region (90 • -100 • E, 0 • S−10 • N) (Saji et al., 1999). For ENSO we adopt the Niño-3.4 index, defined as the area-averaged SSTAs over 170 • W−120 • W, 5 • S−5 • N. To examine the strength of oceanic linkages between the Pacific and Indian oceans and the relationship between ENSO and IOD, we perform a correlation analysis using the depth anomalies in the 20 • C isotherm (Z20), a proxy for the tropical thermocline depth, in the ITF region (area average over 120 • E−131.5 • E, 7.5 • -8.5 • S) (Bracco et al., 2005) or the sea level difference between the western Pacific (WP; 0-10 • N; 125 • E−145 • E) and the Eastern Indian Ocean (EIO; 10 • S−20 • S; 110 • E−130 • E) (Mayer et al., 2018).

THE ENTROPY QUANTIFIER: RECURRENCE PLOTS AND INFORMATION ENTROPY
Given a dynamical system, several measures of complexity have been proposed to distinguish regular (e.g., periodic), chaotic and random behaviors in the time-series that describes it. Among those complexity quantifiers, the most common are Lyapunov exponents, fractal dimensions, scaling exponents, divergence rates and entropies (Manneville, 1990;Shi, 2007). Most of these quantifiers work very well in the case of low dimensional dynamical systems but their application is complicated for noisy time series coming from high-dimensional real-world systems. Bandt and Pompe (2002) proposed the permutation entropy quantifier to address this issue, and along with minor modifications this tool has been used in climate science for the analysis of a proxy record of ENSO spanning the Holocene (Saco et al., 2010) and for ice core records by Garland et al. (2019). The permutation entropy depends on two parameters and at least for low dimensional systems (for example a Logistic Map) has a good, but not optimal, correlation with the Lyapunov exponents of the system under investigation.
Here we adopt instead an entropy quantifier based on the distribution of microstates in a recurrence plot (Prado et al., 2020). All complexity quantifiers are based on fundamental phase space properties of ergodic dynamical systems, such as trajectory recurrence (Poincaré, 1890;Cvitanović et al., 2016) and this recurrence can be visualized, for a given time-series, using the recurrence plot (RP) method introduced by Eckmann et al. (1987). Given a trajectory x i in a dynamical system in its ddimensional phase space at time i, its RP is given by an N × N matrix of 1 and 0 such as: where ε is the threshold distance and defines the neighborhood of a state x i , is the Heaviside function, · is a norm (Euclidean distance in our case) and N is the number of states considered. The analysis of structures (such as diagonal, vertical or horizontal lines) in a RP is known as recurrence quantification analysis (RQA) and has found numerous applications in biology, neuroscience, physics, geosciences and economics, among other disciplines (Webber and Marwan, 2015). In essence, given a time series it is possible to reconstruct a state space representation using the embedding theorem (Takens, 1981), and then compute the recurrence plot (which is a matrix with "1" if two states are recurrent, "0" if they are not recurrent). Most methods to estimate complexity of time series are directly related to RQA (Marwan et al., 2007), recurrence network analysis  and information theory (Bandt and Pompe, 2002;Balasis et al., 2013). Some of these quantifiers have been applied to climate science to analyze regime shifts and tipping points in single time-series (e.g., Donges et al., 2011). All these methods require "phase space reconstruction, " as mentioned, or, in other words, they require embedding the time series in a higher dimensional space. Phase space reconstruction is computationally time consuming and sensitive to both the choice of parameters and to the system's noise and dimensionality (Gilpin, 2020), and is therefore not advisable for real-world spatiotemporal fields as it may lead to spurious results (Marwan, 2011;Riedl, 2013).
To remedy these issues, Corso et al. (2018) proposed a recurrence entropy quantifier that does not require phase space reconstruction and can be safely applied to fields of any dimension and complexity. The idea behind is that the entropy of a time series can be computed by the probability of occurrence of microstates in its RP. A microstate of size M is defined as an M × M matrix inside the RP. The total number of configurations of 1 and 0 in a microstate of size M is M * = 2 M 2 and is possible to define a probability of occurrence P k of a microstate k as P k = n k M * , with n k being the number of occurrences of the k-th microstate in the RP. The information entropy of the time series is then given by: It is possible to normalize the recurrence entropy by its maximum possible value M 2 ln 2 which corresponds to the case in which all microstates appear with the same probability P k = 1 M * ; in this case S = 0 (1) implies perfect predictability (unpredictability) of the system dynamics.
The number of admissible microstates M * grows exponentially as a function of M but Corso et al. (2018) showed that only few of them populate the RP. It is therefore reasonable to consider only a random subsample M of the possible microstates in the RP and expect a rapid convergence to S (M * ) (see Supplementary Material). In this case, the entropy is normalized by ln M.
In layman's terms, the information entropy based on the recurrence of microstates can be computed in two steps by calculating firstly the recurrence plot and then the Shannon Entropy based on the distribution of microstates. If the Euclidean distance between states x(t = a) and x(t = b) is less than epsilon then they are considered neighbors and described by "1" in the recurrence plot matrix, vice versa they are described by "0." The recurrence entropy, in whichever way it is calculated, depends on the choice of the distance threshold ε. Several heuristics have been proposed for selecting ε : ∼5% of the maximal phase space diameter (Mindlin and Gilmore, 1992); no more than 10% of the mean (or maximum) phase space diameter (Koebbe and Mayer-Kress, 1992;Zbilut and Webber, 1992); a value ensuring a recurrence point density of ∼1% (Zbilut et al., 2002). More recently, Prado et al. (2020) relied on the "maximum entropy principle" to free the analysis from having to select ε. Given a RP and microstates of size M, they examined the dependence of the recurrence entropy S on ε and found that, typically, the curve S(ε) has a well-defined maximum, S Max , that does not vary significantly for a range of ε values. In simple deterministic and stochastic systems, the maximum entropy S Max calculated by Prado and colleagues is highly correlated with the Lyapunov exponent, further improving the complexity quantifier proposed by Corso et al. (2018). S Max is therefore our preferred choice. Technical details for the heuristic used to determine it and an example of the convergence when considering only a subsample M of all possible microstates in the RP can be found in the Supplementary Material.
The information entropy so computed measures the complexity of the time series examined, and the two terms, entropy and complexity will be used interchangeably in the reminder of the paper. The higher the entropy, the more complex and less predictable is the system, and vice versa. Examples from simple systems such as a logistic map (May, 1976) can be found in both Corso et al. (2018) and Prado et al. (2020). In the Supplementary Material we quantify the information entropy of the Lorenz system (Lorenz, 1963), while a first application to noisy, high-dimensional data can be found in Falasca et al. (2020) in which the authors used the entropy quantifier to identify abrupt regime shifts in paleoclimate simulations.
The information entropy when applied to a climate field quantifies the degree of complexity of a given region (or grid point) and depends only on the size of a microstate M. In the next section we will explore its sensitivity from M = 2 to M = 8. In summary, for a given climate field X(t), in this work we compute the information entropy and its variability for all the time series in the spatial grid to derive the spatiotemporal entropy field S X (t) without having to rely on embedding.

ENTROPY IN REANALYSIS DATA, CMIP5 HISTORICAL RUNS AND CMIP5 PROJECTIONS
The entropy fields for the SST in the SODA/ERA-Interim reanalysis and two representative CMIP5 models, CCSM4 and HadGEM2-ES are shown in Figure 2 using monthly average fields and all months, and for different values of M. With this first analysis we explore the robustness of metric to the size of the microstates. A microstate of size M = 3, for example, quantifies the recurrence between x(t 1 ), x(t 1 +1 month), x(t 1 +2 months) and another random sequence x(t 2 ), x(t 2 +1 month), x(t 2 +2 months), while if M = 2 the recurrence will be evaluated between x(t 1 ), x(t 1 +1 month) and another, shorter, random sequence x(t 2 ), x(t 2 +1 month), and so on.
The patterns in the entropy plots do not vary significantly with M, but the entropy values do, with generally lower predictability (higher entropy values) and smaller chances of recurrence for increasing M, as to be expected given that we are evaluating the recurrence of a longer (in time) microstate, and very high, unstructured, predictability if M is too small. This dependence can be exploited for studies focusing on different time-scales and output frequency, for example by using daily data for evaluating subseasonal predictability.
When all months are considered, in SODA/ERA-Interim the highest predictability is found in the central tropical Pacific in the ENSO impacted area, as expected. Upwelling regions, both along South America in the cold tongue and in the Arabian Sea, are characterized by higher complexity near the coasts that in offshore waters, and the IOD and SIOD action regions appear slightly more predictable than the remaining IO.
CCSM4 overestimates predictability (underestimates entropy) nearly everywhere. In terms of patterns, it is close to the observed ones, but it does not capture the lower than surrounding predictability of the cold tongue in the Pacific and in the western Arabian Sea upwelling, and underestimates slightly that of the ocean region to the west of the tip of India, that participate in the IOD dynamics. HadGEM2-ES, on the other hand, does not capture the ENSO predictability potential around the Equator (10 • N−10 • S), and overestimates predictability in the upwelling systems in the Pacific and to the west of Australia.
The usefulness of the entropy is also in highlighting changes in predictability among months, seasons and subseasonal time scales if daily or higher frequency data are available. Given the monthly averages used in this study, we show next the entropy in boreal spring (March to May, MAM) (Figure 3), when complexity is expected to be higher due to the spring predictability barrier to ENSO, and in the extended fall season, August to November (ASON) (Figure 4), when the impacts on the rainfall variability in the continents surrounding the Indian Ocean are greatest and the variance explained by the IOD is largest.
In both figures the entropy is shown in SODA/ERA-Interim and all models using M = 3, chosen in light of the seasonal scope of our analysis (3 or 4 instead of 12 months in each year). In all CMIP5 models the complexity is higher in boreal spring in the Pacific, indicating that they generally capture the spring barrier to ENSO predictability (e.g., Torrence and Webster, 1998;Duan and Wei, 2013). Little agreement is found among models on the spring predictability of IO.
In boreal fall, models with a realistically shaped-but underestimated-entropy in the ENSO region, especially in the southern hemisphere (CCSM4, CanESM2, MPI) tend to overestimate the predictability in the IO, west of Sumatra for CCSM4 and west of Australia in CanESM2 and MPI. This bias was not as evident in spring. Other models share the bias in the IO but also have high (generally too high) predictability in the tropical Pacific and lower than observed at the equator (ACCESS1.0 and 1.3, CMCC-CESM). CNRM-CM5 reproduces best the entropy patterns in both basins. GFDL-CM3 underestimates the complexity of the equatorial Pacific, while the IPSL model, in both its version, is characterized by an ENSO pattern protruding too far west into the warm pool area and very high predictability north and south of the equator in the western Pacific. The remaining models tend to underestimate the SODA/ERA predictability in the Indo-Pacific. In the IO, many models underestimate the entropy in the 15 • S−35 • S band and overestimate it along the equatorial ocean, missing the IODrelated predictability and displaying little intra-model agreement on the overall patterns.
Since the IOD discovery, Saji et al. (1999) reported the existence of a correlation between the IOD and ENSO indices. This correlation reaches 0.62 in the ASON season over the period considered for the reanalysis, consistent with the SODA/ERA-Interim entropy pattern, which suggests some predictability potential in the eastern, equatorial and part of the western IO. In the Arabian Sea, on the other hand, the complexity tends to be larger, likely due to the energetic mesoscale field characterizing this upwelling system in fall.
Since correlation does not imply causation, many studies have questioned whether the IOD can occur independently of ENSO and by which mechanism (e.g., Allan et al., 2001;Ashok et al., 2003;Gualdi et al., 2003;Saji and Yamagata, 2003). The general consensus is that the IOD is an ocean-atmospheric coupled mode of climate variability intrinsic to the IO that can be excited by ENSO (see also Webster et al., 1999;. As a result, its predictability does depend on ENSO through an atmospheric teleconnection and possibly through an oceanic bridge. The atmospheric teleconnection is indisputable and is achieved through changes in the Walker circulation over the western Pacific Ocean during the development of ENSO events (e.g., Lau and Nath, 2003;Annamalai et al., 2005;Lau et al., 2005;Behera et al., 2006;Kug et al., 2006;Izumo et al., 2010;Kajtar et al., 2015). Some studies suggest that the ENSO-IOD interaction may also be modulated by an oceanic bridge through changes in the Indonesian Throughflow (ITF) transport between the two basins (Bracco et al., 2005;Yuan et al., 2011;Zhou et al., 2015). The ENSO modulations of the thermocline depth in the Pacific Ocean, and particularly in the Warm Pool region, propagate through the ITF to the northwestern Australia coast and then into the IO (Cai et al., 2005;Behera et al., 2006). The ITF signal is then transported to the southeastern tropical IO (SETIO) by coastal Kelvin waves that develop off south Java (Sprintall et al., 1999). Modeling studies found indeed that closing the ITF annihilates the ENSO-IOD relationship (Wajsowicz and Schneider, 2001;Bracco et al., 2005;Song et al., 2007;Santoso et al., 2011;Kajtar et al., 2015) and van Sebille et al. (2014) confirmed that the ITF transport correlates with ENSO using an eddy-permitting ocean model. In section Entropy and the ENSO-IOD Relation we will briefly investigate these atmospheric and oceanic connections in light of the large discrepancies found among models and between models and reanalysis in the entropy field.
We focus next on future projections. The multi-model mean consensus on the evolution of ENSO under global warming in CMIP5 is toward a strengthening of its amplitude and permanent El-Niño-like conditions (IPCC, 2014;Cai et al., 2015). At the same time the IO warms up but not uniformly. The predicted warming differs significantly among models, but in most the western Indian Ocean warms more than the east side (Di Nezio et al., 2020). In all models but the two versions of GISS-E2, CNRM-CM5 and INMCM4, the warming pattern over the IO resembles that of a positive IOD, which is associated with above normal precipitations over East Africa in fall and greater chances of fires in Australia and Indonesia (not shown). Despite agreement among many models in the warming patterns, no robust behavior is found in the model representation of tropical SST predictability in the RCP8.5 global warming scenarios (Figure 5). In the Pacific, ACCESS1.0, CanESM, IPSL-CM5A-LR are characterized by a decrease in entropy, the opposite is verified in ACCESS1.3, CMCC-CESM, GFDL-CM3, both versions of the GISS model, HadGEM2-ES, IPSL-CM5A-MR, MPI-ESM-LR and INMCM4. The remaining models are almost invariant. In the Indian Ocean the behavior may differ from the Pacific. The entropy decreases in the ACCESS and IPSL runs, all four of them, in CanESM, MRI and in the GFDL model south of the Equator, increases in the CNRM-CM5 run, which had the closest representation to the reanalysis in the historical period, and in the GISS-E2-R model, and is nearly unchanged in the rest of the cases.
The historical and future entropies suggest that many CMIP5 models misrepresent both ENSO and the IOD in boreal fall, and their relationship, and generally underestimate the complexity of the Indian Ocean SST variability, which is neither a mere response to ENSO or independent of it.

ENTROPY AND THE ENSO-IOD RELATION
In this section, we first discuss how the entropy measure relates to more traditional analysis methods and then we delve into the model biases in the ENSO-IOD representation. We focus on the historical period, for which reanalysis products are available for comparison.
The information entropy is a non-linear time series analysis tool that allow for extracting and accounting for non-linear information that cannot be resolved using traditional linear methods, such as Empirical Orthogonal Function (EOFs), or power spectra. At the same time, the entropy bias can be anticipated in part from a traditional EOF analysis of SST anomalies (Figure 6). Indeed, the entropy identifies major problems in the representation of climate variability patterns, as done by the EOFs, if these biases influence recurrence. Additionally, the information entropy introduces information on the complexity of the climate modes in different regions. Too large complexity in the equatorial Pacific, for example, is indicative of a structural problem common to several models which is not (or not only) linked to the modeled SST ENSO patterns per se, but to the representation of equatorial dynamics both in the atmosphere and in the ocean. The Indo-Pacific EOF patterns are generally more in agreement with reanalysis data than the entropy field. ENSO is often modeled with the observed strength and maximum load at the equator even when it extends too far into the warm pool region (e.g., ACCESS1.3, CCSM4, GFDL, GISS-E2-H, INMCM4, both IPSL versions and MPI) and has too much strength in the central and western portion of the Pacific basin compared to the eastern one (Bellenger et al., 2014;Chen et al., 2017). It is however too strong in CMCC-CESM and CanESM, and too weak in HadGEM2, MRI and INMCM4. In the reanalysis, the SST anomaly patterns in the IO are characterized by a modest negative signal into the SETIO region that extends into the IO from the warm pool region, and by a positive signal elsewhere with slightly larger intensity in the Arabian Sea. The SETIO pattern is often underestimated (this is the case for CanESM and HadGEM-ES, both versions of the GISS and IPSL models), and is stronger than observed but equatorially confined in GFDL. The positive signal, on the other hand, is overestimated in CMCC, GFLD, IPSL and MPI. The EOF patterns, again, do not correspond to the entropy ones, especially south of 15 • S, where the SST predictability is often greatly overestimated by the models.
The complexity of ENSO is also-but not exclusively-linked to its power spectrum and therefore its periodic recurrence in time. Figure 7 compares the spectra among models and the reanalysis based on the (monthly) Nino3.4 index. The bias in the modeled spectra is reflected in that of the entropy field along the equatorial Pacific, with the caveat that several models overestimate the complexity of ENSO at the equator but display highly predictable (and therefore too regular) dynamics away from it, especially in the southern hemisphere (see e.g., the ACCESS models). Noticeably, models with too much power that spreads across multiple time scales maintain relatively low predictability, being the last related to the repetitiveness of the microstates (see e.g., CMCC-CESM). Given the interval considered, low entropy characterizes also models with ENSO power on time scales that are too long (e.g., HadGEM2-ES, for which ENSO peaks at 8-10 years).
In Figure 8 we show the variance of the Niño-3.4 and IOD indices across the four seasons. In many models, the variances of the two indices are comparable in magnitude, while being characterized by an approximately 3:1 ratio in the observational dataset. Additionally, the Niño-3.4 -IOD crosscorrelation (Supplementary Figure 3), which is linked to the directionality of the connections between the two modes and therefore the predictability potential across the basins, indicates that the Niño-3.4 variability lags, or co-occurs with, the IOD one in most models, while it leads by 2 months in the reanalysis. CNRM-CM5 emerges again as the model with the closest representation to the reanalysis. Cross-correlations between the Nino3.4 index and the SETIO and WTIO indices separately point to the eastern region of the IO as source of the bias (not shown).
To better highlight what limits the reliability of model climate predictions in the IO, we next examine briefly atmospheric and oceanic connections induced by developing ENSO events affect the IO and vice versa. We use standard correlation and regression analysis and our investigation does not aim at being exhaustive, but simply identifies possible biases that are shared among many models and appear to influence how predictability is simulated.

Atmospheric Teleconnection
The sea level pressure anomalies (SLPAs) over the SETIO region can be used as a proxy for the strength of atmospheric changes to the Walker circulation induced by ENSO over this area. A simple cross-correlation analysis shows that the modeled SLPa-ENSO relationship is captured by all models (Figure 9). The same is verified for the WTIO region (not shown). Interestingly CNRM-CM5 underestimates the signal, which is better captured by the GISS runs, CanESM or GFDL.
From this simple analysis, we can infer that the bias in representing the IOD-ENSO interaction is not caused, at a first order, by a deficient representation of the atmospheric teleconnection from the Pacific into the IO. The IOD can independently modify the Walker circulation through air-sea interactions in SETIO (Izumo et al., 2010) but the ENSO influence on the Walker circulation is captured relatively well in all models. We will show later that the same cannot be said if we zoom onto the ITF region.

Oceanic Bridge
ENSO influences the IO ocean circulation by modulating the transport of warm and fresh water from the Pacific Warm Pool region into the IO via the ITF. Positive transport anomalies occur during La Niña events and vice versa in El Niño years (Meyers, 1996). England and Huang (2005) found the ITF transport, defined as the depth integrated transport over the whole water column at 8 • S, 120 • -131.5 • E, to be anticorrelated with the Niño-3 index (SSTa averaged over 150 • W−90 • W, 5 • S−5 • N), with a maximum of −0.32 when the ITF lags ENSO by 9 months. Bracco et al. (2005) focused instead on the variability of the 20 • C isotherm (Z20) that may more directly impact the upwelling in the SETIO region and found that that the fall correlation  between the Z20 anomalies and Niño-3.4 reached −0.62 over the period 1958-2002. In SODA, the maximum (minimum) transport occurs in boreal spring (fall), and the thermocline depth oscillates between 123 and 135 m. In the models, the annual cycle varies greatly in both amplitude and phase (Supplementary Figure 4). Several models have a thermocline that is 30 m or more too deep. Many, including INMCM4, GISS-E2-R, GISS-E2-H, MPI and CANESM, display an out of phase seasonal cycle (minimum in spring and maximum in fall) or shaped with two peaks. CMCC-CESM stands out for having a thermocline that is too shallow (about 20 m shallower than in the reanalysis in fall). CCSM, GFDL and IPSL-MR provide the closest representation to that in the reanalysis, with CCSM being the most realistic.
The cross-correlations between the Z20 anomalies averaged over the ITF region and Nino3.4 indices are shown in Figure 10. In the reanalysis, the maximum anticorrelation (−0.71) is found for the Niño-3.4 index lagging the ITF thermocline by 1 month. Statistically significant anticorrelations are found in all models but INMCM. However, the modeled anticorrelations are consistently maximized at about 3-5 months lags. The relationship between the Niño-3.4 index and the ITF is misrepresented by the models also when using the sea level difference between the western Pacific (WP; 0-10 • N; 125 • E−145 • E) and the Eastern Indian Ocean (EIO; 10 • S−20 • S; 110 • E−130 • E) as ITF proxy (Mayer et al., 2018). Correlation coefficients are smaller than for Z20 but significant in the reanalyses, while the linkage is missed entirely in most models (Supplementary Figure 5). This is indicative of problems in the representation of both oceanic and atmospheric processes at the equator.
Finally, the regression maps of the SST over the tropics and of subsurface temperatures averages between 5 • N and 5 • S onto the IOD index are shown in Figures 11, 12. At the surface the link between the IOD and the ENSO anomalies are captured by most models but INMCM4 and underestimated by ACCESS1.3, HadGEM2-ES, followed by CCSM4, MPI and MRI and CNRM-CM5. At the ocean surbsurface, on the other hand, the regression patterns are generally strong underneath the warm pool, where they have the same sign of the SETIO region, but underestimated in the central and eastern Pacific. At the subsurface CNRM-CM5 emerges as closest to the reanalysis in the pattern representation. Modeled IODs somewhat mimic the characteristics of the unseasonal IOD events identified in the observations by Du et al. (2013), but for their seasonality. The unseasonal IOD are relevant in late spring and summer since the 1970's, and resemble the CMPI5 ones by being mostly independent of ENSO and forced by (too strong for fall in the CMIP5 case) equatorial winds. This figure also points to the differences in the representation of the ITF bathymetry among the models.
The representation of the relative seasonal evolution of the two modes identified in the cross-correlations and the greater independence of the modeled IOD from the Pacific evolution emerge as problematic in all models. Our brief analysis extends the results by  that highlighted other systematic biases in the IO and in particular (a) the western IO thermocline being deeper than the eastern IO while the opposite is verified in the observations so that the wave dynamics that transport the IOD temperature anomalies are not reproduced correctly; (b) slightly warmer SST in the western IO than eastern IO, while a much warmer eastern IO is observed; (c) stronger than observed climatological easterly winds over the equatorial Pacific.

SUMMARY
Weather prediction is reliable to <10 days due to the atmosphere chaotic behavior, the sensitivity on the initial conditions and model limitations. Climate prediction, on the other hand, has the potential to extend to longer time scales due to the modulations exerted by its slowly varying component, the ocean. Climate modes that manifest as non-linear oscillations, such as the El Niño-Southern Oscillation, the Pacific Decadal Oscillation, the Atlantic Multi-Decadal Oscillation, the Indian Ocean Dipole, etc. further contribute to the extended range of potential predictability of the climate system.
Global warming is posed to modify these modes of natural variability that control much of the predictability at seasonal to interannual scales. Examples include the intensification of the IOD (Di Nezio et al., 2020), and a higher frequency of Central Pacific El Niños (Freund et al., 2019) observed in the last decade and projected to continue in the future, and shifts in its global-reaching teleconnections following the expansion of the Hadley cells (Kang and Lu, 2012). In light of the societal impacts that these modes exert by modulating surface temperatures, storm frequency, droughts and floods, their current and future predictability potential is of the uttermost importance. It remains, however, elusive due to the complexity of the climate system, the biases in climate models, the challenges implicit in quantifying it.
In this work, we introduced a new method, building upon tools developed within the non-linear dynamical systems community, to quantify predictability in terms of information entropy. The information entropy of a climate field quantifies the degree of complexity of a given region or grid point in terms of recurrence (how many times a phase space trajectory visits roughly the same area in the phase space over a given period) and evaluate how it varies over time. Given that the evolution of a climate field such as SST results from interactions of many fields, from atmospheric winds and heat fluxes to ocean currents, its space phase evolution accounts for all these interactions, linear and non-linear.
Most studies employing a comparable non-linear data analysis technique in climate science have focused on a limited number of 1-dimensional time series (i.e., paleoclimate proxies-one or few at the time-or their modeled equivalent) (Donges et al., , 2015. The possibility to use it on large, multi-dimensional fields has been opened by Corso et al. (2018) by introducing a new entropy quantifier applicable to d-dimensional fields using the recurrent plot (RP) technique (Eckmann et al., 1987) without the need for embedding or phase space reconstruction. We applied it to the SST field in the tropical Pacific and Indian Ocean, comparing CMIP5 models and reanalysis data, over 39 years in the historical period and in the last 30 years of the XXI century RCP8.5 projections. We focused on the boreal fall season, when the ENSO signature is strong and the IOD has maximum variance in the Indian Ocean. Understanding the ENSO-IOD relation and the influence of the IOD on all countries facing the Indian Ocean is essential to improve prediction of extreme rainfall and droughts over Australia, India, Indonesia and East Africa. In these highly populated regions downscaled regional experiments (Hashizume et al., 2012) are limited by the representation of the large scale atmospheric and oceanic circulation.
Our analysis highlights the limitations that coupled climate models still face in addressing climate predictability and its potential changes. No robust signal is found in the models, in both basins, as to how predictability will evolve in the future when considering the projections for the last 30 years of the XXI century. In the historical period, very few models capture both pattern and intensity of the entropy signal of the reanalysis, and CNRM-CM5 outperforms all others in this regard. Interestingly, in an intercomparison performed using complex network analysis (Fountalis et al., 2015), CNRM-CM5 outperformed all other CMIP5 models also in the representation of rainfall patterns and their connectivities, but performed as well as several others in representing the (global) SST network. This again points to the power of the entropy quantifier as a measure of topological similarity linked to predictive power. The entropy quantifier, in summary, is a computationally efficient performance metric for coupled climate models, capable of capturing dynamical characteristics, both linear and non-linear, that goes beyond the climatology and local variability of the field under investigation, and delves into its topological properties.
We also showed that a biased representation of coupled equatorial dynamics and of the atmospheric and subsurface oceanic bridge between the Pacific and Indian Oceans via the ITF contributes to the poor representation of the Indo-Pacific entropy in fall. Incidentally, we verified that the intra-model variability is relatively small in the historical runs in the statistics presented, from the entropy to the cross-correlations, and each model behavior is consistent across ensemble members. Overall, many models struggle at the equator, in both basins, and display unrealistic regular dynamics in the IO, especially south of the equator and/or in the Arabian Sea.
Our results exemplify how information entropy may contribute a new powerful tool to investigate the potential predictability of the climate system. More work is needed to explore its relevance on time scales other than interannual, using higher or lower frequency time series as input to explore, for example, subseasonal phenomena o changes across millennia, and the entropy usefulness in analyzing fields other than SSTs. In the near future, we aim at applying it across time scales and with high frequency (at least daily) data to quantify how and where climate predictability emerges from the weather noise.