Assessing causal dependencies in climatic indices

We evaluate causal dependencies between thirteen indices that represent large-scale climate patterns (El Nino/Southern Oscillation, the North Atlantic Oscillation, the Pacific Decadal Oscillation, etc.) using a recently proposed approach based on a linear approximation of the transfer entropy. We demonstrate that this methodology identifies causal relations that are well-known, as well as it uncovers some relations which, to the best of our knowledge, have not yet been reported in the literature. We also identify significant changes in causal dependencies that have occurred in the last three decades.


Introduction
Identifying significant causal dependencies directly from time-series analysis is a challenging problem with practical applications in all fields of science and technology. Main challenges are to distinguish causality from correlation, and to distinguish direct from indirect causality (considering, for example, three processes X, Y and Z, X can directly affect Z, or it can indirectly affect Z because X affects Y and in turn, Y affects Z). In climate research, an additional challenge is to identify non-trivial causalities, i.e., those that are not due to geographical proximity nor due to common seasonal effects. Many approaches for causal inference have been proposed in the literature, based of model reconstruction, information theory, phase-space reconstruction, nonlinear symbolic analysis, machine learning, etc. (Granger 1969;Schreiber 2000;Baccala and Sameshima 2001;Eichler et al. 2003;Sugihara 2012;Runge 2018;Siyang Leng et al. 2020;Subramaniyam et al. 2021;Huang et al. 2020). These approaches have different mathematical hypotheses and computational complexity, and they differ considerably in their capability to detect genuine couplings (Krakovska et al. 2018). Nevertheless, they have been extensively used in climate research and we mention here just a few examples (for a comprehensive review we refer the reader to Runge (2019)).
In Mosedale et al. (2006) Granger causality (GC) (Granger 1969) was used to study the effect of sea surface temperatures (SSTs) on the North Atlantic Oscillation (NAO) on daily time scale. It was found that the SST tripole index provides additional predictive information for the NAO than that available by using only past values of NAO, i.e., it was found that the SST tripole is Granger causal for the NAO.
In Liang (2014) a formula for evaluating the rate of information flow from one series to another was proposed and applied to study the relation between El Niño and the Indian Ocean Dipole (IOD). It was found that the causality is asymmetric: El Niño tends to stabilize IOD, while the effect of IOD on El Niño is more uncertain. The rate of information flow was recently used in Docquier et al. (2022) to analyze the cause-effect relationships between Arctic sea ice and its potential drivers. It was found that changes in Arctic sea ice are mainly driven by air and sea-surface temperatures and ocean heat transport. The rate of information flow was also used in Vannitsem and Liang (2022) to assess directional dependencies of seven climate indices and three local Belgium time series, on time scales from a month to five years. Among other results, the authors found that the Arctic Oscillation plays a key role on the dynamics of the North Pacific (NP) and 1 3 North Atlantic (NA) on short (monthly) time scale, while NAO plays a role on long time scales (several years). These findings are consistent with earlier work (Vannitsem and Ekelmans 2018) using convergent cross mapping (CCM) (Sugihara 2012), as it was found that the Niño3.4 region influences the NA dynamics through its annual climatological cycle, the atmosphere over the NP region forces the NA region in a monthly basis and there is mutual influence of NP and NA on longer (interannual) timescales.
CCM was also used in Zhang et al. (2019) to examine causality between the Northern Hemisphere annular mode (NAM) and the wintertime surface air temperature (SAT) over Northeast Asia. It was found that NAM information is encoded in SAT data but not the other way around, indicating a causal link in the direction NAM→SAT.
In Tirabassi et al. (2015) GC was used to study the airsea interaction in the South Atlantic Convergence Zone (SACZ) and it was shown that GC allows distinguishing different regimes: there are years in which the forcing is mainly directed from the atmosphere to the ocean, years in which the ocean forces the atmosphere, years in which the influence is mutual and years in which the coupling is not significant.
In Tirabassi et al. (2017) renormalised partial directed coherence (that is a form of Granger causality in the frequency domain (Baccala and Sameshima 2001)) and directed partial correlation (that is an alternative approach to quantify GC in the time domain (Eichler et al. 2003)) were used to analyze the NINO3.4-Monsoon interaction and also, the air-sea interaction in the SACZ. While a clear NINO3.4-Monsoon mutual interaction was found, regarding SACZ, the results were not as clear but reinforced the evidence that, in some years, SACZ may be forced by SST.
In Nowack et al. (2020) an approach based on iteratively testing for conditional independent relationships among time series (Runge 2018(Runge , 2019) was used to analyze sea level pressure data from a large set of climate model simulations and, as a proxy for observations, meteorological reanalyses. It was shown that the causal networks obtained offer an objective pathway for model evaluation and model inter-comparison. In Di Capua et al. (2020) this conditional independence-based approach was used to assess the direction of causal interactions and to quantify the relative causal effects of tropical, mid-latitude and internal drivers on the Indian summer monsoon rainfall variability.
In Manshour et al. (2021) the transfer entropy, TE (Schreiber 2000), which is a form of conditional mutual information, CMI (Hlavackova-Schindler et al. 2007), was used to study causality and information transfer between the solar wind and the magnetosphere-ionosphere system. The causal relations uncovered could be described in terms of linear time-delayed information transfer, with delays of 10 to 30 minutes.
In Deza et al. (2015) a directionality index defined in terms of the transfer entropies TE Y→X and TE X→Y , that allows to measure the net flow of information between two processes X and Y, was used as an estimator of the net direction of information flow in a global dataset of surface air temperature anomalies. The structures found recovered some well-known climate variability patterns, such as atmospheric waves in the extratropics and longer range events in the tropics.
Recently, two of us used a measure, referred to as pseudo transfer entropy (pTE) (Silini and Masoller 2021) that is a linear estimator of the TE, as it is the analytical expression of TE for Gaussian processes. Previous works have used the Gaussian TE expression to study, using the wavelet transform, causality across rainfall time scales (Molini et al. 2010), and to study, using the Hilbert transform, phase-amplitude coupling in a century long record of data of daily surface air temperature from various European locations (Palus 2014).
The pTE measure has two parameters that have to be properly selected: the order of the model used to represent the data, and the lag time of information transfer. In Silini and Masoller (2021) the pTE results were validated using model generated data (where the underlying causality was known) and then, it was applied to the analysis of the tele-connection between the Central Pacific and the Indian Ocean, as inferred from the bivariate analysis of the monthly-sampled time series of NINO3.4 and All India Rainfall (AIR) indices. While pTE and GC only detected the dominant causality (NINO3.4 → AIR), TE detected both.
Keeping in mind that the pTE methodology might only return strong causal links (direct or indirect, while it might miss weak or nonlinear but significant causalities), in this work we aim to extend the previous study and analyze a set of indices that represent main patterns of climate variability. Our goal is to detect significant causalities regardless of whether they are due to direct or indirect teleconnections. Our motivation is that, in either case, the gained knowledge can contribute to improve the forecast of a climatic pattern, such as ENSO, by taking into account information of the subset of indices that are causally related to the index that represents the pattern. We are also interested in shedding light on how causal interactions might have changed over the last decades.
This paper is organized as follows. Section 2 describes the datasets analyzed; Sect. 3 describes the methodology used; Sect. 4 presents the results and Sect. 5 presents the discussion of the results.

Data
We focus this study in the climatic indices described below. They are monthly sampled timeseries and are freely accessible at the NOAA website (https:// psl. noaa. gov/ data/ clima teind ices/ list/), with the exception of the All Indian Rainfall index, which is available via the Indian Institute of Tropical Meteorology (https:// www. tropm et. res. in/). AIR: All Indian Rainfall. The area-weighted integral of the rainfall measured by the Indian national network of rain gauges.
AMO: Atlantic Multidecadal Oscillation. The detrended area-weighted average over the North Atlantic, from the equator up to 70N, of the sea surface temperature (SST) anomalies from the Kaplan SST dataset (Kaplan et al. 1998;Reynolds and Smith 1994).
GMT: The Global Mean Temperature anomaly as computed by NASA/GISS. The anomaly is computed with respect to the period 1951-1980.
HURR : The total number of hurricanes or named tropical storms in a given month in the Atlantic region.
NAO: The North Atlantic Oscillation. The north-south dipole of pressure anomalies over the North Atlantic, with one center over Greenland and the other center of opposite sign between 35N and 40N. SOI: Southern Oscillation Index. The standardized difference in surface air pressure between Tahiti and Darwin. The SOI is a proxy of the strength of the Walker circulation and it's strictly related to ENSO.
The various indices span different regions and focus on different variables. variability of the ocean and the atmosphere on different spatio-temporal scales, with particular attention to the tropical belt.
The majority of the timeseries span six decades, overlapping in the period 1950-2016 (i.e., in 792 data points). In this period, the timeseries are depicted in Fig. 1.
The various indices display different spectral properties. Several have a defined periodic component, either seasonal, as in the case of rainfall and storm indices, or longer, like the case of the QBO. Some display trends (GMT), others slow non-linear oscillations (NINO34). The high frequencies as well are very heterogeneous.
All these complex spectral properties influence the indices' distributions. Generally, we observe skewed distributions and hints of multimodality.
AIR, HURR, Sahel and NP indices display a strong seasonal component (Fig. 1). To avoid spurious signals, we removed the seasonality subtracting the mean of every different month from these four indices. The other indices have already been constructed after removing the seasonal cycle of the variables used.
Since the indices under consideration span a broad range of values, we standardized their distribution by removing a linear trend and rescaling the series to have unitary variance. Even if information-related quantities such as pTE, transfer entropy, or mutual information in principle do not depend on the variables gauges, the rescaling to unitary variance is considered good practice for their numerical computation. The resulting post-processed time series and their values distributions are depicted in Fig. 2.

pseudo Transfer Entropy
We analyze the causality structure of the set of indices using bivariate analysis, i.e., we address the problem of determining whether one index influences another without considering the possible influence of a third index that could mediate the interaction.
To measure the bivariate causality between two indices, we use the pseudo transfer entropy (pTE) (Silini and Masoller 2021), a method that has recently proven to be a computationally fast alternative to traditional transfer entropy (TE) (Schreiber 2000). pTE is a linear estimator of TE, as it is the analytical expression of TE for Gaussian processes and therefore, it is equivalent to Granger causality, GC (Barnett et al. 2009). Analogously to TE, the pTE measures the transfer of information from Y at time t to X at time t + conditional to the information flowing from X at time t to X at time t + , providing a mean to assess causal relationships between two processes.
The pTE from series Y to X is calculated as (Silini and Masoller 2021) where n+1 is the vector containing the future elements of the time series X, k n is the matrix containing the k past values of X, ⨁ stands for the horizontal concatenation and Σ[A] is the covariance of matrix A.
We note that pTE is an asymmetric measure, thus in general pTE Y→X ≠ pTE X→Y . We also note that if the past of Y does not influence the future of X, pTE Y→X = 0 , else pTE Y→X > 0 . To determine whether a small pTE Y→X value (1) is consistent with the null hypothesis that Y has no influence on the future of X, a significant test needs to be performed, as described in the next section. The need of a significance test is a limitation of pTE analysis, as the results can depend on the type of surrogate data used to perform the test and on the significance criterion used (Silini and Masoller 2021). On the other hand this is a well-known limitation, shared by other methods used in the literature. Another limitation is the fact that pTE is a linear estimator and thus, it can miss some nonlinear causalities. Additionally, when applied to non-Gaussian data, pTE may return fake causalities. On the other hand, a main advantage is that Eq.(1) allows an efficient calculation in comparison with the original TE formulation (Schreiber 2000). This is important when a large number of links need to be tested. A detailed comparison was performed in Silini and Masoller (2021), and here we mention just as an example that for the analysis of two time series of 500 data points each, The embedding dimension, k, has to be selected before carrying on the calculation. There are different possibilities to determine its optimal value. Here we model X as an autoregressive process and fix the model order, k, minimizing the Bayesian information criterion (BIC) score.
In Table 1, we report the values of k that minimize the BIC score for each index. As expected, for the majority of the atmospheric indices, a short-memory AR process is sufficient. In particular, we can see that the order of NAO is 1, and this is consistent with the fact that a first order approximation of NAO is red noise. Ocean indices are overall modelled with a higher order AR process, indicating a longer memory. Nevertheless, the PDO and TSA indices are modelled as AR(1) processes, even though they are dominated by long-term variations, indicating that in these cases  autocorrelation is mediated by the effect of short lags. The high k = 8 value for NTA can be explained in part by considering that, when generating the index, the SST anomalies were smoothed by a 3-months running mean.

Statistical significance analysis
Once the pTE between two indices is computed, we have to address whether it is significant or not. Unlike the case of cross-correlation, the null model of X being independent of the past of Y doesn't allow the analytical calculation of p-values for the pTE distribution. For this reason, we have to rely on surrogate analysis to understand if a pTE value is significantly different from zero. Surrogate timeseries can be obtained from real-world data through different kinds of manipulation (Lancaster et al. 2018). Surrogate time series should retain all the properties of the original timeseries with the exception of the one we are interested in. In the case of causal relationships, we want surrogates that are independent from each other while preserving the autocorrelation function of the original time series. In fact, preserving the autocorrelation ensures that we preserve the dependence of the timeseries on itself. To achieve this we employed an algorithm known as iterative amplitude adjusted Fourier transform (IAAFT) (Schreiber andSchmitz 1996, 2000), which preserves both the amplitude distribution and the power spectrum of the original series.
From the original dataset, we generated N = 1000 independent surrogate datasets using IAAFT. From the surrogate datasets we obtain N surrogate measures of pTE between each pair of timeseries. Thus, the quantiles of the surrogate pTEs can be viewed as significance threshold for the measured pTE values. In the following, we considered a pTE value significant if it falls within the highest 1% of its surrogates distribution.

Long-term causality variation
Measuring the variation of pTE between the first and second half of the dataset allows us to explore possible long-term variation in the index interactions.
We first divided the time series of each index in two nonoverlapping segments with the same length, corresponding to the periods pre-1984 and post-1984. All the time series are monthly-sampled and have 804 data points (that cover the period 1950-2016), therefore, we have 402 data points in each segment. Then, we measured the pTE between each pair of indices in the two halves and we computed the difference.
To address the significance of such differences we rely again on surrogate analysis. In particular, we created 1000 surrogate datasets for each one of the two halves and we computed the pTE for each pair of indices in the 2000 surrogates. For each pTE of the second half of surrogates, we sampled 100 pTEs from the first half and we computed the average difference. Hence, we end up with 1000 surrogate average differences for each pair. If the original difference is within the highest or lowest 1% of the surrogates, depending on its sign, than we consider that difference significant.

Results
In Fig. 3a we report an example of pTE calculation focusing on NINO34 as "forcing" node, considering significant pTEs for different values of . Most of the results match with the current knowledge regarding ENSO dynamics. We can observe a 5-months cutoff in the NINO34 → SOI interaction, which is in line with the ENSO build-up time scale when the ocean and the atmosphere are coupled. Moreover, a maximum of around 4 months in the NINO34 → NTA is expected too, given that their interaction is mediated by heat fluxes: because of its thermal inertia, the SST of the ocean boundary layer changes in a time scale of roughly three months, producing the pTE delayed maximum. From this perspective, the behavior of the AMO is analogous. The AIR and HURR have a 1-month impact, which is reasonable given that the interaction is mediated directly by the atmosphere. It is interesting to note that the HURR index has a small but significant pTE tail up to = 3 , which may result from indirect interactions mediated by the NTA.
Results of some indices are, however, unexpected to a degree. For the PDO and the NP, we would have expected a behavior more similar to the NTA one. Instead, we observe in Fig. 3a relatively high pTEs up to 4 months. We interpret this as due to the fact that the local air-sea interaction increases the persistence of the remotely forced ENSO signal. In contrast, the pTE values for NAO (Fig. 3b) show a very rapid decrease with , indicating interactions on a much shorter time scale.
In the following, for every pair of indices, we calculated the pTE and its significance for = 1, 3 and 6, that is one month, one season and half a year into the future. In Fig. 4, we report the significant connections between the various indices. For better clarity, results are displayed both as a network and an adjacency matrix. In the network representation, directed connections are represented by arrows. In particular, we draw a link only if the value of the pTE between two indices is significant for at least one value of .
We investigate the role of in the pTE values in Fig. 5, again representing connections by directed arrows. The arrows' colors represent the value of for which the pTE is the highest, while the arrows' width represents this maximum value. We can observe that, in general, connections have relatively low . Also, the values of pTE for = 6 seems to be lower than for other lags.
The full network picture allows to visualize the global structure; however, the different links can not be clearly distinguished. For the sake of clarity, we select a subset of key indices (AMO, NTA, PDO, and the NINO34/ SOI pair), and report in Fig. 6 the forward and inward links separately. The pivotal role of ENSO in the climate network is evident, with numerous forward connections tying NINO34 and SOI to the vast majority of the studied indices.
Finally, we report in Fig. 7 the significant variations of pTE values between pre-1984 and post-1984. In Fig. 7, we observe that the variations are heterogeneous, with main changes being an increase in the strength of the NP-PDO link and the decrease of the strength of the PDO-GMT link, and of several links that affect AMO. The discussion of the uncovered links, and their variations, is presented the next section. Fig. 3 Influence of NINO34 a and NAO b to a subset of indices for various lags . We report only those indices for which at least one value of pTE was significant. The vertical axis shows the relative pTE value, which is the relative difference between the measured pTE and the significance threshold determined by using the surrogate analysis (see Sect. 3 for details). Values of 0 mean that the pTE for those lags was not significant (e.g.

Discussion
While some of the uncovered connections are well-known, others are either previously unknown or are possible false positives. In the latter case, a common driver of two indices could produce a significant amount of shared information, inducing an apparent connection. Let's start with the connection between ENSO and the north Pacific indices (NP and PDO). It is well known that ENSO induces an atmospheric teleconnection that modulates the Aleutian low over the North Pacific (Wang et al. 2012). As the Aleutian low is characterized by the NP index, the ENSO → NP causality is well understood.
Changes in the surface winds associated with the Aleutian low in turn induce SST anomalies in the north Pacific thus affecting the PDO (NP → PDO).
On the other hand, it is unlikely that the number of tropical storms in the north Atlantic (HURR) can drive NAO, TSA, ENSO and NP, as suggested by the results. Instead, these causalities likely result as consequence of complex interactions among the different atmospheric and oceanic phenomena characterized by these indices. For example, it is well known that during El Niño the number of hurricanes in the north Atlantic decreases because of enhanced vertical shear (Klotzbach 2011) as found by pTE. At the same time a warm NTA, which is also influenced by ENSO (Chang et al. 2003), favours the development of tropical storms (Pérez-Alarcón et al. 2021). Similar complex interactions explain the links of the AIR index.
The lag is a crucial parameter that has a large impact on the analysis. The global surface air temperature warms up by about 0.1 • C during an El Niño event, with a lag of about 6 months (Trenberth et al. 2002). This causality is detected in the analysis, although ENSO → GMT is maximum at lag 1. At longer lags we find the opposite causality (GMT → SOI at lag 3 and NINO34 at lag 1), which may be understood as consequence of the persistence of the ENSO events that last between 6 and 9 months.
The causality identified QBO→ NAO in our analysis has been reported in the literature to occur during boreal winter (Andrews et al. 2019;Marshall and Scaife 2009), although the link had been assessed as relatively weak. We found the largest causality for a lag of 6 months, suggesting that the mechanism through which the QBO affects the NAO may last more than one season.
It is well known that the NAO is the main driver of the sea surface temperature anomalies over the tropical north Atlantic mainly through changes in surface heat fluxes  (Visbeck et al. 2001). However, we don't find this causality, probably due to the index used to describe the Atlantic SST. As mentioned above the NTA index uses SST anomalies that have been smoothed thus filtering out the response to NAO. An influence of NAO on AMO and NTA was found in Vannitsem and Liang (2022), and a possible reason for the difference with our results is that the NTA index used in Vannitsem and Liang (2022) is different from the one we use here, as our index is smoothed by a 3-month window.
On the other hand, the analysis detected the NTA → NAO connection at lag 1, which is consistent with the literature that shows that SST in the tropical north Atlantic can induce atmospheric teleconnections that project onto NAO (Okumura et al. 2001).
Another index that presents links with several other phenomena is the TSA. As in the case of the tropical north Atlantic, it is likely that the TSA → NAO link is direct, as the SSTa in the south Atlantic can control the position of the Intertropical Convergence Zone which could promote the development of a teleconnection to the north Atlantic (Okumura et al. 2001). Interestingly, connections with other indices occur with lags of 3 or 6 months, suggesting that some of these links are indirect. The tropical south Atlantic is known to influence the equatorial Pacific through changes in the Walker circulation with a lag of several months (Rodríguez-Fonseca et al. 2009), consistent with our results. Thus, we hypothesize that the connections of the TSA with the PDO and HURR indices occurr via the ENSO influence.
Our analysis also shows that the impact of TSA on ENSO has grown in recent decades (see Fig. 7), in agreement with the literature (Rodríguez-Fonseca et al. 2009).
Looking at Fig. 7, the strongest and most consistent signal is the change in causality between SOI (and Nino34) and AMO. That AMO variability influences the ENSO variability is already documented (Levine et al. 2017). While the literature on the link between AMO and ENSO is extensive, to our knowledge, there is no report that this influence is getting stronger, which could have implications for the ENSO predictions. As already mentioned, El Niño warms up the NTA, which is part of the AMO index, therefore it's not surprising that ENSO appears driving AMO. On the other hand, the Atlantic can influence ENSO variability by changing the mean state in the equatorial Pacific, altering the Walker circulation and trades, as pointed out above. From Fig. 7, we infer that the AMO impact on ENSO is increasing while the ENSO impact on AMO is becoming weaker. For longer lags (e.g., = 9 ) only the link AMO → ENSO remains (not shown). The link ENSO → AMO is weak for long leads, while AMO → ENSO can still be strong, due to the different time scales of the phenomena. Figure 7 also shows an increase in the NP → PDO link in the last decades, which may be related to the fact that the ENSO teleconnection to the north Pacific has also increased (NINO34 → NP, SOI → NP). On the other hand the link PDO → NP does not seem to have changed. Combined these results suggest that the SST anomalies in the north Pacific Fig. 7 Significant differences in the causality networks between two time windows, represented both as network and adjacency matrix. The first time window contains the years 1950-1983, while the second correspond to the window 1984-2016. Green (red) links correspond to an increment (reduction) in the causality from the first to the second window. For each connection, we report the largest significant difference across the three studied values of . The value of tau for which the difference is the largest is reported in adjacency matrix. The width and color of the links are proportional to the relative increments (reductions), with respect to the significance threshold found with IAAFT surrogates (p-value = 0.01) have become more dependent on the equatorial Pacific conditions compared to local air-sea interactions.
Regarding the time lags of significant causalities, while NAO is causal for only a month or two, ENSO's causal effects spread up to 6 months. This result could be due to the different memories of the indices since ENSO has a much longer persistence than NAO. IAAFT surrogates preserve the power spectrum of the original series, retaining the autocorrelation structure and mitigating this problem. Nevertheless, it is possible that the different indices memories are still impacting the results, and further surrogate testing is needed to clarify the origin of the link. We conclude the discussion by remarking that, while in our study we have found lags that are, in general, consistent with known lagged physical influences, this may not be the case when dealing with chaotic or nonlinearly coupled systems, and the detection of coupling delay is a challenging problem that has not yet been solved in general terms (Coufal et al. 2017).

Conclusions
We have used the pseudo transfer entropy (pTE), which is a simplified expression of the transfer entropy, to evaluate causal dependencies between thirteen indices that represent large-scale climate pattern. Taken together, our results have unveiled the well-known complexity of the network of interactions and feedback loops, and their interdecadal variations. The majority of the links recovered by our analysis have been documented in the literature and can be explained through known physical mechanisms; however, we have also found undocumented or likely spurious interactions. While it is important for advancing the understanding of our climate to identify the links that represent genuine connections, from a practical standpoint, to improve the forecast of an index variability, the pTE analysis yields useful knowledge because it tells us which signals contain information relevant for the future of another signal. In this way, the pTE represents a useful tool of time series analysis, to identify features that can potentially improve the forecast of the evolution of a climate index. As an example, we have found that the link between AMO and ENSO is becoming stronger, which may be important for ENSO predictability.
Another application of the pTE algorithm is for performing model inter-comparisons, e.g., for contrasting the causal links found in model data with those found in observed data, in order to determine the skill of different climate models in representing the interactions and lags in our climate.
Funding Open Access funding provided thanks to the CRUE-CSIC agreement with Springer Nature. This work received funding from the European Union's Horizon 2020 research and innovation program under Grant Agreement No 813844. C.M. also acknowledges funding by the ICREA ACADEMIA program of Generalitat de Catalunya and Ministerio de Ciencia e Innovacion, Spain, project PID2021-123994NB-C21.

Data availability
The datasets analysed during the current study are freely available at the links https:// psl. noaa. gov/ data/ clima teind ices/ list/ and https:// www. tropm et. res. in/. The code to compute the pTE and the IAAFT surrogates is freely available at https:// github. com/ ricca rdosi lini/ pTE

Conflict of interest
The authors have no relevant financial or non-financial interests to disclose.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.