The Effect of Time Series Distance Functions on Functional Climate Networks

TL;DR

This paper systematically evaluates 29 time series distance functions to understand their impact on the structure of functional climate networks derived from global temperature data, revealing diverse long-distance connection patterns.

Contribution

It provides a comprehensive analysis of how different distance functions influence climate network structures, offering new tools for climate variability analysis.

Findings

Different distance functions produce distinct network structures.

Some functions reveal unique long-distance connection patterns.

Alternative measures can enhance understanding of climate dynamics.

Abstract

Complex network theory provides an important tool for the analysis of complex systems such as the Earth's climate. In this context, functional climate networks can be constructed using a spatiotemporal climate dataset and a suitable time series distance function. The resulting coarse-grained view on climate variability consists of representing distinct areas on the globe (i.e., grid cells) by nodes and connecting pairs of nodes that present similar time series. One fundamental concern when constructing such a functional climate network is the definition of a metric that captures the mutual similarity between time series. Here we study systematically the effect of 29 time series distance functions on functional climate network construction based on global temperature data. We observe that the distance functions previously used in the literature commonly generate very similar networks while alternative ones result in rather distinct network structures and reveal different long-distance connection patterns. These patterns are highly important for the study of climate dynamics since they generally represent pathways for the long-distance transportation of energy and can be used to forecast climate variability on subseasonal to interannual or even decadal scales. Therefore, we propose the measures studied here as alternatives for the analysis of climate variability and to further exploit their complementary capability of capturing different aspects of the underlying dynamics that may help gaining a more holistic empirical understanding of the global climate system.

Figures

Tables

Table 1. Table 1: Time series distance functions used in the experiments

	Distance	PDF	Reference
01	Autocorrelation Coefficients (ACF)		[47]
02	Avg($L_1$, $L_{\infty }$) (avgL1LInf)	✗	[12]
03	Bhattacharyya	✗	[53]
04	Compression-based (CDM)		[51]
05	Complexity Invariant (CID)		[38]
06	Correlation		[13]
07	Cross-correlation		[13]
08	Dice	✗	[54]
09	DISSIM		[39]
10	Dynamic Time Warping (DTW)		[37]
11	Euclidean ($L_2$)		[13]
12	Fourier Coefficient (fourierDist)		[46]
13	Gower	✗	[55]
14	Integrated Periodogram (INTPER)		[48]
15	Jaccard	✗	[13]
16	Kulczynski	✗	[13]
17	Lorentzian	✗	[13]
18	Manhattan		[13]
19	Mutual Information (MI)	✗	[44]
20	Maximal Information Coefficient (MIC)	✗	[45]
21	Normalized Compression (NCD)		[50]
22	Partial Autocorrelation Coefficients (PACF)		[47]
23	Periodogram (PER)		[49]
24	Sorensen	✗	[56]
25	Squared-Euclidean (sqdEuclidean)	✗	[13]
26	Short Time Series (STS)		[40]
27	Tanimoto	✗	[57]
28	Temporal Correlation and Raw Values (CORT)		[43]
29	Wave Hedges	✗	[12]

Table 2. Table 2: Time series distance functions used in this work

Distance	Equation	Ref.
Avg($L_1$, $L_{\infty }$) (avgL1LInf)	$d_{ALL}(P,Q)=\left({\sum }_i^t\left|P_i-Q_i\right|+\max \left|P_i-Q_i\right|\right)/2$	[12]
Bhattacharyya	$d_{BH}(P,Q)=-\ln {\sum }_i^t\sqrt{P_i{Q}_i}$	[53]
Dice	$d_{DC}(P,Q)=\frac{{\sum }_i^t{\left(P_i-Q_i\right)}^2}{{\sum }_i^t{P}_i^2+{\sum }_i^t{Q}_i^2}$	[54]
Gower	$d_{GW}(P,Q)=\frac{1}{t}\sum _i^t\left|P_i-Q_i\right|$	[55]
Jaccard	$d_{JC}(P,Q)=\frac{{\sum }_i^t{\left(P_i-Q_i\right)}^2}{{\sum }_i^t{P}_i^2+{\sum }_i^t{Q}_i^2-{\sum }_i^t{P}_i{Q}_i}$	[13]
Kulczynski	$d_{KS}(P,Q)=\frac{{\sum }_i^t\left|P_i-Q_i\right|}{{\sum }_i^t\min (P_i,Q_i)}$	[13]
Lorentzian	$d_{LO}(P,Q)=\sum _i^t\ln \left(1+\left|P_i-Q_i\right|\right)$	[13]
\pbox30cmMutual Informationa (MI)	$d_{MI}(X,Y)=H\left(X\right)+H\left(Y\right)-2MI(X,Y)$	[44]
	$H\left(X\right)=-{\sum }_i^t{P}_i\log P_i$
	$MI(X,Y)=\sum _i^t\sum _j^{t}p(X_i,Y_j)\log \frac{p(X_i,Y_j)}{P_i{Q}_j}$
\pbox22cm Maximal Information
Coefficientb (MIC)	$d_{MIC}(X,Y)=1-MIC(X,Y)$	[45]
	$M_{xy}=\frac{\max \left(MI\left(G_{xy}\right)\right)}{\log \left(\min (x,y)\right)}$
Sørensen	$d_{SO}(P,Q)=\frac{{\sum }_i^t\left|P_i-Q_i\right|}{{\sum }_i^t\left(P_i+Q_i\right)}$	[56]
\pbox20cmSquared Euclidean
(sqdEuclidean)	$d_{SED}(P,Q)=\sum _i^t{\left(P_i-Q_i\right)}^2$	[13]
Tanimoto	$d_{TM}(P,Q)=\frac{{\sum }_i^t\left(\max (P_i,Q_i)-\min (P_i,Q_i)\right)}{{\sum }_i^t\max (P_i,Q_i)}$	[57]
Wave Hedges	$d_{WH}(P,Q)=\sum _i^t\frac{\left|P_i-Q_i\right|}{\max (P_i,Q_i)}$	[12]

Table 3. Table 3: Time series distance functions used in this work

Distance	Equation	Ref.
\pbox20cmEuclidean	$d_{\text{ED}}(X,Y)=\sqrt{{\sum }_{t=1}^T{\left(X_t-Y_t\right)}^2}$	[13]
Manhattan	$d_{\text{MH}}(X,Y)={\sum }_{t=1}^T\left|X_t-Y_t\right|$	[13]
\pbox22cmComplexity
Invariant
(CID)	$d_{\text{CID}}(X,Y)=cf(X,Y)\cdot d_{\text{ED}}(X,Y)$	[38]
	$cf(X,Y)=\frac{\max (ce\left(X\right),ce\left(Y\right))}{\min (ce\left(X\right),ce\left(Y\right))}$
	$ce\left(X\right)=\sqrt{{\sum }_{t=1}^{T-1}{\left(X_t-X_{t+1}\right)}^2}$
\pbox22cmDynamic
time
warpinga
(DTW)	$d_{\text{DTW}}(X,Y)=dtw\left(i=T,j=T\right)$	[37]
Correlationb	$d_{\text{COR}}(X,Y)=1-\left|{\rho }_{XY}\right|$
	${\rho }_{XY}=\frac{{\sum }_{t=1}^T\left(X_t-\overline{X}\right)\left(Y_t-\overline{Y}\right)}{\sqrt{{\sum }_{t=1}^T{\left(X_t-\overline{X}\right)}^2}\cdot \sqrt{{\sum }_{t=1}^T{\left(Y_t-\overline{Y}\right)}^2}}$	[13]
Cross-correlationb
(crossCor)	$d_{\text{CCOR}}(X,Y)=1-\left(\underset{\tau \in [-{\tau }_{max},{\tau }_{max}]}{\max }\left|{\rho }_{X{Y}_{\tau }}\right|\right)$	[13]
	${\rho }_{X{Y}_{\tau }}=\frac{{\sum }_{t=1}^T\left(X_t-\overline{X}\right)\left(Y_{t-\tau }-\overline{Y}\right)}{\sqrt{{\sum }_{t=1}^T{\left(X_t-\overline{X}\right)}^2}\cdot \sqrt{{\sum }_{t=1}^T{\left(Y_{t-\tau }-\overline{Y}\right)}^2}}$
\pbox22cmAutocorrelation
Coefficientsb
(ACF)
Partial
Autocorrelation
Coefficients
(PACF)	$d_{\text{ACF}}(X,Y)=\sqrt{\sum _{\tau =-{\tau }_{max}}^{{\tau }_{max}}{\left({\rho }_{X{X}_{\tau }}-{\rho }_{Y{Y}_{\tau }}\right)}^2}$	[47]

Table 4. Table 4: Time series distance functions used in this work

Distance	Equation	Ref.
\pbox22cmTemporal Correlation and Raw Values d (CORT)	$d_{\text{CORT}}(X,Y)=\Phi \left(cort(X,Y)\right)\cdot d_{\text{ED}}(X,Y)$	[43]
	$\Phi \left(u\right)=\frac{2}{1+e^{ku}}$
	$cort(X,Y)=\frac{{\sum }_{i=1}^{T-1}\left(X_{t+1}-X_t\right)\left(Y_{t+1}-Y_t\right)}{\sqrt{{\sum }_{i=1}^{T-1}{\left(X_{t+1}-X_t\right)}^2}\sqrt{{\sum }_{i=1}^{T-1}{\left(Y_{t+1}-Y_t\right)}^2}}$
\pbox22cmShort Time
Seriesc (STS)	$d_{\text{STS}}(X,Y)=\sqrt{{\sum }_{t=1}^{T-1}{\left(\left(Y_{t+1}-Y_t\right)-\left(X_{t+1}-X_t\right)\right)}^2}$	[40]
DISSIM	$d_{\text{DISSIM}}(X,Y)=\sum _{t=1}^{T-1}{\int }_t^{t+1}\left|X_t-Y_t\right|dt$	[39]
\pbox22cmNormalized
Compressione
(NCD)	$d_{\text{NCD}}(X,Y)=\frac{C\left(XY\right)-\min (C\left(X\right),C\left(Y\right))}{\max (C\left(X\right),C\left(Y\right))}$	[50]
\pbox22cmCompression-
basedd
(CDM)	$d_{CDM}(X,Y)=\frac{C(X,Y)}{C\left(X\right)+C\left(Y\right)}$	[51]
\pbox22cmPeriodogram
(PER)	$d_{\text{PER}}(X,Y)=\sqrt{\sum _{k=1}^{T/2}{\left(I_X\left({\lambda }_k\right)-I_Y\left({\lambda }_k\right)\right)}^2}$	[49]
	${\lambda }_k=\frac{2\pi k}{T}$
	$I_X\left({\lambda }_k\right)=\frac{1}{T}{\left|{\sum }_{t=1}^T{X}_t{e}^{-i{\lambda }_{k}t}\right|}^2$
\pbox22cmIntegrated
Periodogram
(INTPER)	$d_{\text{INTPER}}(X,Y)={\int }_{-\pi }^{\pi }\left|F_X\left(\lambda \right)-F_Y\left(\lambda \right)\right|d\lambda$	[48]
	$F_X\left({\lambda }_j\right)=\frac{{\sum }_{i=1}^j{I}_X\left({\lambda }_i\right)}{{\sum }_i{I}_X\left({\lambda }_i\right)}$
\pbox22cmFourier
Coefficientf
(fourierDist)	$d_{FO}(X,Y)=\sqrt{\sum _{j=1}^n{\left|x_j-y_j\right|}^2}$	[46]
	$x_j=\sum _{k=1}^T{X}_k\cdot e^{2\pi i\left(k-1\right)\left(j-1\right)/T}$

Table 5. Table 5: Global measures from networks created with percentile p = 0.001 𝑝 0.001 p=0.001

distance	$LC$	$L$	$C$	SWI	$r$	$Q$
acf	25.00 %	39.48	0.93	1.82	0.97	0.78
avgL1LInf	67.39 %	7.77	0.87	9.73	0.95	0.72
bhattacharyya	68.79 %	6.80	0.84	10.98	0.92	0.71
cdm	62.85 %	3.50	0.61	19.71	0.22	0.55
cid	9.32 %	18.01	0.90	3.32	0.94	0.68
correlation	63.54 %	105.45	0.92	0.66	0.95	0.76
cort	21.63 %	53.40	0.90	1.12	0.93	0.67
crossCorr	64.00 %	105.45	0.92	0.66	0.95	0.76
dice	69.00 %	9.15	0.88	7.93	0.95	0.73
dissim	37.27 %	29.64	0.87	2.30	0.93	0.75
dtw	7.84 %	9.85	0.89	6.54	0.94	0.74
euclidean	17.89 %	49.56	0.90	1.20	0.93	0.66
fourierDist	8.45 %	16.38	0.90	3.61	0.92	0.66
gower	65.84 %	8.13	0.89	9.00	0.96	0.73
intper	28.30 %	35.15	0.95	1.96	0.98	0.81
jaccard	69.00 %	9.15	0.88	7.97	0.95	0.73
kulczynski	65.84 %	8.13	0.89	9.03	0.96	0.73
lorentzian	69.05 %	7.97	0.87	9.24	0.95	0.72
manhattan	9.08 %	20.09	0.89	2.94	0.92	0.66
mi	55.31 %	125.36	0.90	0.53	0.95	0.75
mic	69.82 %	70.53	0.94	1.05	0.98	0.81
ncd	62.54 %	3.51	0.61	19.75	0.24	0.55
pacf	34.04 %	43.74	0.91	1.55	0.96	0.74
per	8.00 %	15.26	0.88	3.71	0.90	0.57
sorensen	65.84 %	8.13	0.89	9.02	0.96	0.73
sqdEuclidean	68.85 %	7.41	0.81	10.55	0.90	0.70
sts	41.16 %	63.69	0.91	1.07	0.96	0.77
tanimoto	65.30 %	8.17	0.90	9.02	0.96	0.73
waveHedges	69.05 %	7.97	0.87	9.26	0.95	0.72

Table 6. Table 6: Global measures from networks created with percentile p = 0.01 𝑝 0.01 p=0.01

distance	$LC$	$L$	$C$	SWI	$r$	$Q$
acf	100.00 %	6.10	0.56	8.06	0.73	0.48
avgL1LInf	97.83 %	5.03	0.31	6.29	0.19	0.46
bhattacharyya	96.39 %	4.55	0.31	6.51	0.17	0.44
cdm	97.75 %	2.48	0.11	4.36	-0.14	0.24
cid	100.00 %	6.45	0.76	6.90	0.87	0.66
correlation	100.00 %	10.76	0.76	4.40	0.88	0.64
cort	100.00 %	11.48	0.77	3.97	0.91	0.63
crossCorr	100.00 %	10.74	0.76	4.41	0.88	0.64
dice	98.00 %	4.92	0.35	7.13	0.20	0.47
dissim	99.97 %	4.72	0.48	8.27	0.69	0.67
dtw	95.14 %	7.31	0.62	4.37	0.39	0.56
euclidean	100.00 %	7.11	0.78	6.26	0.88	0.64
fourierDist	100.00 %	6.16	0.77	7.27	0.87	0.71
gower	98.39 %	5.07	0.31	6.66	0.19	0.51
intper	99.49 %	11.17	0.56	4.91	0.75	0.57
jaccard	98.00 %	4.92	0.35	7.13	0.20	0.47
kulczynski	98.39 %	5.07	0.31	6.66	0.19	0.51
lorentzian	98.55 %	5.01	0.31	6.45	0.19	0.46
manhattan	100.00 %	6.58	0.77	6.71	0.88	0.64
mi	100.00 %	10.23	0.76	4.44	0.89	0.65
mic	100.00 %	9.61	0.75	5.01	0.86	0.63
ncd	97.62 %	2.49	0.11	4.37	-0.13	0.24
pacf	99.94 %	4.12	0.65	11.08	0.81	0.59
per	94.00 %	9.40	0.68	2.81	0.51	0.32
sorensen	98.39 %	5.07	0.31	6.66	0.19	0.51
sqdEuclidean	94.34 %	4.80	0.32	6.18	0.20	0.47
sts	99.99 %	8.43	0.73	3.90	0.59	0.58
tanimoto	98.35 %	5.09	0.31	6.71	0.19	0.51
waveHedges	98.55 %	5.01	0.31	6.45	0.19	0.46

Table 7. Table 7: Global measures from networks created with percentile p = 0.05 𝑝 0.05 p=0.05

distance	$LC$	$L$	$C$	SWI	$r$	$Q$
acf	100.00 %	4.17	0.56	3.53	0.59	0.43
avgL1LInf	99.82 %	3.72	0.46	3.36	0.32	0.39
bhattacharyya	99.77 %	3.47	0.46	3.51	0.30	0.40
cdm	99.94 %	2.03	0.21	2.65	-0.16	0.18
cid	100.00 %	2.64	0.53	5.03	0.46	0.40
correlation	100.00 %	3.06	0.65	5.32	0.45	0.56
cort	100.00 %	2.50	0.61	6.86	0.64	0.52
crossCorr	100.00 %	3.02	0.62	5.00	0.42	0.55
dice	100.00 %	3.58	0.47	3.58	0.30	0.41
dissim	100.00 %	2.54	0.44	4.34	0.25	0.32
dtw	99.48 %	3.33	0.58	3.52	0.09	0.33
euclidean	100.00 %	2.47	0.51	5.20	0.25	0.42
fourierDist	100.00 %	2.50	0.49	5.04	0.23	0.40
gower	99.88 %	3.66	0.46	3.47	0.33	0.40
intper	99.95 %	6.74	0.60	2.66	0.64	0.55
jaccard	100.00 %	3.58	0.47	3.58	0.30	0.41
kulczynski	99.88 %	3.66	0.46	3.47	0.33	0.40
lorentzian	99.89 %	3.63	0.46	3.45	0.33	0.43
manhattan	100.00 %	2.53	0.52	5.19	0.28	0.41
mi	100.00 %	2.26	0.53	6.17	0.41	0.53
mic	100.00 %	2.84	0.63	5.38	0.40	0.56
ncd	99.93 %	2.04	0.21	2.65	-0.16	0.18
pacf	99.99 %	2.67	0.41	4.32	0.17	0.43
per	99.00 %	4.58	0.76	2.80	0.26	0.12
sorensen	99.88 %	3.66	0.46	3.47	0.33	0.40
sqdEuclidean	98.49 %	3.81	0.48	3.22	0.31	0.39
sts	100.00 %	3.37	0.69	3.73	0.18	0.27
tanimoto	99.88 %	3.66	0.46	3.47	0.33	0.40
waveHedges	99.89 %	3.63	0.46	3.45	0.33	0.43

Table 8. Table 8: Global measures from networks created with percentile p = 0.1 𝑝 0.1 p=0.1

distance	$LC$	$L$	$C$	SWI	$r$	$Q$
acf	100.00 %	3.44	0.64	2.71	0.61	0.35
avgL1LInf	99.92 %	3.17	0.57	2.73	0.43	0.37
bhattacharyya	99.97 %	2.98	0.56	2.83	0.40	0.36
cdm	99.99 %	1.93	0.29	2.20	-0.16	0.16
cid	100.00 %	2.27	0.53	3.41	0.32	0.29
correlation	100.00 %	2.31	0.58	3.64	0.29	0.41
cort	100.00 %	2.00	0.46	3.67	0.23	0.33
crossCorr	100.00 %	2.37	0.59	3.42	0.24	0.30
dice	100.00 %	3.08	0.57	2.80	0.41	0.39
dissim	100.00 %	2.17	0.48	3.16	0.08	0.28
dtw	99.88 %	2.58	0.61	3.02	-0.01	0.23
euclidean	100.00 %	2.09	0.49	3.37	0.04	0.31
fourierDist	100.00 %	2.16	0.49	3.32	0.05	0.30
gower	99.90 %	3.14	0.56	2.80	0.43	0.37
intper	100.00 %	5.12	0.67	2.12	0.70	0.48
jaccard	100.00 %	3.08	0.57	2.80	0.41	0.39
kulczynski	99.90 %	3.14	0.56	2.80	0.43	0.37
lorentzian	99.90 %	3.11	0.57	2.77	0.44	0.37
manhattan	100.00 %	2.13	0.51	3.41	0.07	0.33
mi	100.00 %	1.90	0.42	3.45	0.26	0.35
mic	100.00 %	2.07	0.55	3.81	0.22	0.38
ncd	99.99 %	1.94	0.30	2.20	-0.16	0.16
pacf	100.00 %	2.35	0.46	3.06	0.20	0.33
per	100.00 %	3.41	0.80	2.71	0.24	0.09
sorensen	99.90 %	3.14	0.56	2.80	0.43	0.37
sqdEuclidean	99.41 %	3.31	0.58	2.57	0.42	0.36
sts	100.00 %	2.56	0.62	2.93	-0.17	0.19
tanimoto	99.90 %	3.14	0.56	2.80	0.43	0.37
waveHedges	99.90 %	3.11	0.57	2.77	0.44	0.37

Equations

$$d_{L_{p}}(X,Y)=\Bigg{(}\sum_{i=1}^{t}(X_{i}-Y_{i})^{p}\Bigg{)}^{\frac{1}{p}}.$$

$$\mathbf{A}=\Theta (\tau -\mathbf{D}),$$

$$d_H(A^G,A^{G^{\prime }})=\sum _{i,j}\frac{1}{2}\mathrm{xor}(A_{i,j}^G,A_{i,j}^{G^{\prime }}),$$

$$d_{Haversine}((\phi_{1},\lambda_{1}),(\phi_{2},\lambda_{2}))=2r\arcsin\Bigg{(}\sqrt{\sin^{2}\bigg{(}\frac{\phi_{2}-\phi_{1}}{2}\bigg{)}+\cos(\phi_{1})\cos(\phi_{2})\sin^{2}\bigg{(}\frac{\lambda_{2}-\lambda_{1}}{2}\bigg{)}}\Bigg{)},$$

$$SWI=\frac{C(G)/C(G_{rand})}{L(G)/L(G_{rand})},$$

Full text

1 Introduction

Statistical climatology provides a great variety of tools for studying the dynamics of meteorological variables across time and space. As a recent approach complementing more traditional analysis techniques, during the last decades complex network tools have emerged as an alternative to statistically analyze and model complex systems such as the Earth’s climate [1, 2, 3, 4, 5, 6, 7, 8]. In this context, the resulting networks are called functional climate networks and can be constructed in different ways. The most common procedure consists of dividing a spatiotemporal climate dataset into grid cells. In this case, network nodes represent the cells of the grid and links are established if the time series of climate variability associated with different cells fulfil a certain criterion of statistical similarity. The main goal of such a climate network description is to highlight the essential interrelationships among the spatially distributed grid points that govern certain relevant aspects of climate variability regionally or even globally.

One question that has not been systematically addressed so far in the context of constructing functional climate networks concerns the choice of a metric that captures similarities between two time series. The majority of previous works has used linear statistical association measures like covariance or Pearson correlation coefficient with or without time-lag [4, 5, 7, 9, 10]. With these measures, it has been possible to construct climate networks and detect interesting climatological patterns, while they have been conceptually restricted to statistical similarities captured by the concept of correlation or straightforward nonlinear generalizations thereof like mutual information [4, 11]. However, there exist a plethora of other distance functions that could be equally applied to the construction of climate networks [12, 13, 14, 15]. For simplicity reasons, we use here the term “distance function” to generally refer to these measures even though some of them are no distance functions in the strict mathematical sense. Some of these functions can capture statistical similarities in climatological variables that correlations cannot detect. Thus, the central question here is whether those alternative distance functions can generate climate networks that highlight different climatological phenomena than the “classical” ones and consequently reveal other meaningful patterns.

In this paper, we evaluate how different choices of time series distance functions affect the properties of the resulting functional climate networks. We are particularly interested in the topological differences among the resulting networks and in understanding the meaning of the climatological patterns captured by them. We use 29 different distance functions and apply complex network theory [16, 17, 18] to analyze the networks constructed based on global near-surface air temperature data. In the construction process, we employ a parameter $p$ that controls the number of edges (i.e., the edge density) for all the networks. This procedure generates networks with the same number of edges. When the edge density is low ($p\to 0$), the resulting networks have only a few edges connecting the most similar time series, while those networks tend to exhibit disconnected components. Conversely, high edge densities ($p\to 1$) generate densely connected networks, but many of the resulting edges connect pairs of time series that are not quite similar. Choosing an intermediate value of $p$ is hence desirable, since it preserves only those edges that represent the most similar time series and avoids the creation of strongly fragmented networks [19].

We start by studying the influence of the edge density $p$ on the networks. For the considered network size (i.e., number of grid points), we observe a transition phase in the interval $0.0004\leq p\leq 0.1$ from disconnected components to connected networks. In this interval, the resulting networks connect only the most similar time series while the networks are still not highly fragmented. In this range of edge densities, the resulting climate networks exhibit a (likely geometrically induced [20, 21, 22, 23]) small-world effect along with degree distributions that decay faster than a power law. Regarding the network topology, we find that the correlation, cross-correlation, mutual information and maximal information coefficient generate similar networks. In contrast, alternative distance functions generate very different network structures. We also observe variations in the connectivity among geographically distant nodes (teleconnections) which we compare with those previously observed in the literature. These long-distance edges, like such emanating from the El Niño–Southern Oscillation (ENSO) region in the eastern tropical Pacific, are important because they represent the transport of energy and atmospheric wave propagation on global scales [24, 25, 26]. Such teleconnections are also practically relevant due to their capability to anticipate climatological events and enhanced forecast skill from subseasonal up to even decadal time scales [27, 28]. Therefore, our results are significant for climate sciences since the identified teleconnection patterns capture distinct types of underlying dynamics. We here show that alternative distance functions permit the study of climate variability from complementary perspectives and could thereby contribute to a better understanding of the climate system. We present these results in details in the following sections.

We organize the remainder of this paper as follows. In Section 2, we introduce the methodological framework employed in this paper. We summarize the relevant concepts of functional climate network construction and time series distance functions, along with the data and specific methodology used to obtain our results. In Section 3, the results of our study are presented in three parts focusing on global network features, network similarities and teleconnection patterns, respectively. Finally, in Section 4, we present some conclusions drawn from our reported findings.

2 Materials and methods

2.1 Functional climate networks

Complx network theory has been applied for studying many different types of complex systems. In the context of the Earth’s climate system, networks can be used to represent the essential information provided by a spatiotemporal gridded climate dataset. In functional climate networks, every grid cell from a given dataset corresponds to a node while only the most similar pairs of nodes are connected in the network (Fig. 1). This approach raises two natural questions: How to measure similarity between time series and which threshold value for this measure to use for identifying two time series as similar and, hence, connected? Both choices have a direct effect on the resulting network topology. The majority of the existing literature [4, 7, 9, 10, 29] has used the absolute value of the Pearson correlation coefficient (PCC), either as instantaneous or lagged correlation, thereby integrating information on both, strong positive and strong negative correlations.

The first application of network science in climate studies was presented by Tsonis et al. [7]. In that work, the authors used wind flow data from the NCEP/NCAR reanalysis dataset [30] and the PCC ($>0.5$) at lag zero to measure similarity between time series. They observed that the resulting network presented an apparent small-world property and was divided into two subgraphs, one of them representing the tropics and the other the extratropical regions. In the tropics, most nodes were densely interconnected, whereas the higher latitudes presented only some distinct highly connected nodes (hubs) and possessed characteristics resembling a scale-free network. The conclusion of the authors at that point was that the subnetwork corresponding to the equatorial region acts as an agent connecting the two hemispheres and thereby allowing for information to travel between them.

Yamasaki et al. [29] used globally distributed temperature time series and observed that the presence of El Niño conditions in the eastern tropical Pacific significantly affects the structure of the resulting network. The authors showed that while El Niño does not necessarily change the temperature in different parts of the world significantly, many network edges are broken under its influence, so that the remaining ones serve as a measure for characterizing this phenomenon. Their tentative conclusion was that the network connections encompass information that had not been captured by previously employed more traditional statistical analyses, which suggests that complex network theory can reveal patterns not observed before by established techniques from statistical climatology. Several later studies have reinforced this idea. For example, Berezin et al. [31] used temperature and geopotential height time series to construct functional climate networks and observed that the relationship between their fluctuations in different geographical regions presents a robust network pattern, while the common belief prior to their study was rather that these variations were not stable and difficult to predict. Steinhaeuser et al. [6] showed that the attributes extracted from communities in climate networks can be used as predictive climatic indices and present statistically better results than traditional clustering methods. Boers et al. [3] demonstrated the possibility to predict climatic extremes using complex networks. Functional climate networks have also been used to investigate large-scale circulation patterns, climate variability modes, their teleconnections, and for model intercomparison [10, 29, 32, 33].

The majority of the aforementioned studies has employed correlation or covariance as distance functions. This implies that these studies have been limited to similarities captured by correlations among the variables under study. In this paper, we show that other distance functions can also be used for the construction of climate networks, an idea that was similarly proposed earlier by Pelan et al. [34]. In the latter work, the authors used six different distance functions and analyzed the effect of edge density on the topology and the clustering predictive capability of the networks. Here, we thoroughly expand this setting by using in total 29 different distance functions and examining the topological differences between the resulting networks.

2.2 Time series distance functions

A time series distance function is a measure of how different two time series are. For the sake of brevity, the term “distance function” will be used here to refer to a time series distance function in the following. Other types of distance functions, e.g., network distance functions, will be referred to with their complete terms in the text. A plethora of distance functions have been proposed in the literature, which can be broadly classified into four main categories [14]: shape-based, edit-based, feature-based and structure-based.

In the following, we provide an overview on some time series distance function concepts that could be used to generate climate networks, without attempting to be complete in this task. Notably, we only consider measures that make use of all time series values, thereby excluding event synchrony measures like event synchronization [3, 28] or event coincidence analysis [35] that have been widely used for studying spatiotemporal patterns of extreme climate events from a complex network perspective. Moreover, we also discard measures of phase synchronization [36] as well as such based on the physical concept of causality and its statistical implementation in terms of coupled stochastic process models or information theory (e.g., Granger causality strength, transfer entropy, and similar measures) from the presented analysis.

Shape-based distance functions.

This kind of measures compare the overall shape of two time series. The most famous measures are $L_{p}$ distances, defined as

[TABLE]

When $p=2$, the $L_{2}$ distance is the well-known Euclidean distance. Other common choices are the $L_{1}$ and $L_{\infty}$ distances that correspond to the Manhattan and Chebyshev distances, respectively. The $L_{p}$ measures are parameter-free and intuitive but fail to capture similarity in many cases, for example, when two time series have similar shapes but are not aligned. This problem is illustrated in Fig 2. $L_{p}$ measures fail in this case because they compare values at fixed positions in the series; they are therefore also called lock-step measures. To solve this conceptual limitation, some alternative elastic measures have been proposed to allow for time warping and provide a more flexible comparison.

One of the most famous concepts for defining elastic measures is Dynamic Time Warping (DTW) [37]. This approach uses a warping path to align two series before their comparison, which is a sequence of adjacency matrix indices that defines a mapping, and the optimal path is the one that minimizes the global warping cost. The Complexity Invariant Distance (CID) [38] is another shape-based measure that involves the Euclidean distance corrected by a complexity estimate of the series. The DISSIM [39] and Short Time Series (STS) [40] distances can deal with time series with different sampling rates by considering linear approximations of the series. DISSIM measures the difference by calculating the integral of the Euclidean distance between them while the STS measure considers the difference of the slopes between measurement points.

Edit-based distance functions.

The idea behind these distance functions is commonly used to compare strings (Levenshtein distance) and consists of counting the minimum number of character insertions, deletions or substitutions to transform one string into another. The Edit Distance for Real Sequences (EDR) [41] considers a threshold value that defines if two real values match or not. The difference is measured by the minimum number of operations to transform one series into another. EDR permits gaps between the two matched subsequences, but assigns penalties according to the lengths of the gaps. The Longest Common Subsequence (LCSS) distance [42] uses an approach similar to the EDR. This measure permits time warping, as DTW, but also allows for gaps in comparison. Different from EDR, the LCSS does not assign penalties to the gaps.

Feature-based distance functions.

These measures extract descriptive features from the time series and compare them. The Pearson correlation coefficient (PCC) $\rho$ is probably the most famous measure of this type. The PCC is defined as the covariance of the two time series divided by the product of their standard deviations. The Temporal Correlation and Raw Values (CORT) measure [43] modulates a conventional distance (like Euclidean or DTW) according to the temporal correlation. Mutual information (MI) is another well known nonlinear measure defined as the information that two time series share [44]. The maximal information coefficient (MIC) considers that if there exits a relationship between two time series, then it is possible to draw a grid on the scatter plot of these two series that encapsulates that relationship [45]. MIC automatically finds the best grid and returns the largest possible MI. Other examples of feature-based measures are based on cross-correlation functions [13], Fourier coefficients [46], auto-correlation coefficients [47] and periodograms [48, 49].

Structure-based distance functions.

Different from feature-based measures, these functions identify and compare high-level structures in the series. Most of them use parametric models to represent time series, e.g., auto-regressive moving average (ARMA) or Hidden Markov Models (HMM). Compression-based measures (CDM) [50, 51] can also be considered to belong to the structure-based distance functions. These measures use compressors like gzip or bzip2 to compress a time series. The distance between two time series is proportional to the difference between the sizes of the compressed files.

2.3 Network generation and analysis

We use the National Center for Environmental Prediction and National Center for Atmospheric Research (NCEP/NCAR) reanalysis dataset [30] of monthly mean near-surface air temperature records on a global regular grid with $2.5\degree\times 2.5\degree$ resolution in latitude and longitude, spanning the period from 1948 to 2016. The resulting data set is composed of 10,512 time series with 816 values each. Prior to our analysis, we remove the seasonal component in each time series by using an additive decomposition by moving averages [52]. We consider here an additive model because the seasonal fluctuations are relatively constant over time. Without removing seasonality, the time series distance functions would show rather trivial patterns reflecting the dominant annual cycle component, while the focus of our present work is on finding similarities in other time series components such as trends or low-frequency variability. For mathematical convenience, we finally normalize all time series to values between 0 and 1.

After preprocessing the data, we calculate the distance between every pair of time series. This step results in a distance matrix $\mathbf{D}$ whose element $D_{ij}$ contains the distance between the series $X_{i}$ and $X_{j}$. For this process, we consider each of the 29 distance functions listed in Tab. 1. The definitions of all distance functions used in this paper can be found in Appendix A in Tab. 2, 3 and A.

For the sake of coherence, we use exclusively distance functions in this paper instead of explicitly considering similarity measures as well. Notably, some functions, like the Pearson correlation coefficient ($\rho$) and mutual information (MI), are actually similarity functions that return high values when two time series are similar. Instead of using $\rho$ directly, we define an associated distance form $d_{\text{COR}}(X,Y)=1-|\rho_{XY}|$. Following the same approach as in previous works [4], we use here the absolute value of the Pearson correlation coefficient because both large negative and large positive correlation values indicate strong statistical interdependence. We apply the same strategy for other similarity measures defined in the range $[0,1]$. The cross-correlation distance is calculated by taking the largest absolute correlation distance for lags in the interval $[-24,+24]$ months. We considered a maximum delay of two years to capture associations delayed by trimesters, semesters, or a year. Conversely, we did not consider longer lags to avoid the drastic reduction in the number of values in the lagged time series. The mutual information (MI) is a similarity function defined in the range $[0,+\infty)$ but commonly bound from above by a maximum value that is specific to the particular estimator used. Instead of using MI itself, we therefore employ here an MI-based distance measure [44] with values in the unit interval $[0,1]$.

Some distance functions measure the dissimilarity by comparing the probability density functions (PDF) of two time series (measures with the PDF column checked in Tab. 1). For these measures, we follow the approach presented in [12] that consists of discretizing every normalized series into 64 bins and dividing the frequencies by the length of the series. The ranges for all the bins are exactly the same because all series were previously normalized to unit range.

After calculating the distance matrix $\mathbf{D}$ for every distance function, we normalize it to simplify the comparison between different measures. Then, we use $\mathbf{D}$ to construct a binary adjacency matrix $\mathbf{A}$ by applying a threshold value $\tau$ as

[TABLE]

where $\Theta(\bullet)$ is the Heaviside step function. Choosing a proper threshold value $\tau$ is fundamental for the network construction process to provide meaningful results. Higher values ($\tau\to 1$) will generate densely connected networks while lower ones ($\tau\to 0$) create fragmented networks. For each distance function, we calculate the distance matrix $\mathbf{D}$ and use the $p^{th}$ distance percentile as the threshold value ($\tau$) to build the network. The distance percentile $p$ is a single parameter that controls the edge density of all networks. Advantages of this approach include the unneccessity to choose a distinct $\tau$ for each distance function and the creation of networks with the same number of edges. Thus, we avoid the topological comparison of networks with different numbers of edges.

As a result of the described process, we obtain a single undirected and unweighted network for every distance function, all of them exhibiting the same number of edges. We compare these networks by using the Hamming distance [58], which measures the number of edges that have to be added or deleted to transform one network into another. Given two labelled networks $G$ and $G^{\prime}$, where $A^{G}$ and $A^{G^{\prime}}$ are their respective adjacency matrices, the Hamming distance is defined as

[TABLE]

where $\mathrm{xor}(\bullet,\bullet)$ is the logical exclusive-or operator.

For defining spatial distances between grid points on a spherical grid, we use the Haversine distance

[TABLE]

where $\phi$ and $\lambda$ correspond to the latitude and longitude of the two locations respectively, and $r$ is the Earth’s radius ($r\approx 6378137$ m).

As a particular scalar characteristic of the obtained networks, we consider the smallworldness index [59] defined as

[TABLE]

where $C(\bullet)$ and $L(\bullet)$ are the respective global clustering coefficient and average shortest path length of the graph $G$ and a correspondent random graph $G_{rand}$, respectively. Considering that this measure uses a probabilistic procedure to generate the random graph benchmark, we repeat this process 50 times and report the mean value taken over all realizations.

One of the concepts used in this paper is a teleconnection, which we define here as the most similar pairs of time series whose grid cells have great-circle distances longer than $5000$ km [60]. By contrast, we disregard short-distance connections since they tend to be similar because of the spatial continuity of the studied temperature field commonly resulting in small time series distances between close nodes. In our analysis, we use the Haversine distance to compute the great-circle distance between two coordinates (Eq. 4). Specifically, we use the center coordinates of each grid cell to define this distance. For keeping the number of results interpretable, we also limit the number of distance functions considered. Instead of studying the networks generated by all the 29 distance functions, we limit ourselves to eight since many distance functions generate similar networks. We have chosen these distance functions by applying hierarchical clustering (Ward’s method) to the Hamming distances between the adjacency matrices of all combinations of networks. After obtaining the eight clusters, we select the centroid network exhibiting the shortest average distance to the other networks in a given cluster.

The choice of $p$ does not change substantially the teleconnections since this is a lower-bound threshold by definition. It means that, for $p=0.01$ or higher values, the strongest links are the same and they all lead to the same teleconnections. Since we show the 500 strongest links (shortest distances), $p$ would just have an impact if it was so small that the network would had less than 500 links. But in this case, $p$ is so small that the network is very disconnected and does not properly represent the system. Therefore, the choice of $p$ is not critical for determining the teleconnections but it has an impact on the networks clustering and on the cluster network centroids. In this paper, we focus on the network centroids as major representatives of the time series distance functions in the cluster but we include the teleconnections for all networks in the Appendix D for the interested reader.

3 Results and discussion

In the following, we evaluate the topological differences of the networks constructed using different time series distance functions. For a review on complex network theory and the full definitions of all measures used throughout this paper, we refer the interested reader to existing review papers of this subject [16, 17, 18].

3.1 Edge density effects

In our experiments, we generate and compare networks with a given edge density $p$. To obtain a better understanding of the effect of this parameter, we first study the corresponding behavior of some important global network measures for different distance functions, which will provide the basis for selecting appropriate values of $p$ for the following analyses. Our results are summarized in Fig. 3.

When the edge density is very low ($p\lesssim 0.0004$), only a few extremely similar pairs of time series are connected, resulting in a very small number of edges and many disconnected components. In technical terms, we have not yet reached the percolation threshold of the network, which can be seen in terms of the very small size of the largest component along with the large number of network components. In this regime, we have a low average path length $L$ (distance between connected nodes), since the majority of nodes belonging to the same component are in fact directly connected. For the same reason, we also find a high global clustering coefficient $C$ (density of closed triangles in the network) and low network modularity $Q$ (heterogeneity associated with community structures). At $p\sim\mathcal{O}(10^{-3})$, we observe a marked change of the network connectivity resulting in the emergence of a giant component along with a marked reduction in the overall number of network components, which is associated with the percolation transition. Notably, this transition occurs at an edge density that is about one order of magnitude larger than for an Erdös-Rényi random graph of the same size ($p_{c}\sim 1/N$), pointing to the fact that edges are not generated at random, but respect the spatial embedding of the network and possibly also tend to cluster in space.

Around the percolation transition, we observe a sharp peak of $L$ and a marked decrease of $C$ at $p\sim\mathcal{O}(10^{-3})$. Both effects are directly related with the increasing connectivity of the network along with the emergence of a giant component. Notably, more and more nodes become (indirectly) connected in the network and thereby contribute to the calculation of both measures. During this phase (about $0.0004\lesssim p\lesssim 0.1$, nodes tend to connect with others of a similar degree (positive assortativity), except for the network created using the CDM distance function that is weakly disassortative (lowest line in Fig. 3F). When $p\approx 0.01$, the networks are almost connected (with only very few isolated nodes) but still sparse. When $p\gtrsim 0.01$, more and more closed triangles appear, and $C$ increases towards its theoretical maximum value ($C=1$) obtained for a fully connected network ($p=1$). On the other hand, the modularity decreases when the network becomes dense since there are no well defined communities anymore. When $p\gtrsim 0.01$, the networks show apparent small-world features (high $C$ and low $L$), which however can be explained by the spatial embedding of the graph as emphasized in several previous studies [20, 21, 22].

In summary, the interval $0.0004\lesssim p\lesssim 0.1$ corresponds to a transition from a disconnected to a strongly connected climate network. After passing this transition phase, the number of edges can further increase until the networks become fully connected ($p=1$). Networks at high edge densities are however not practically interesting since they connect almost all pairs of time series, including those that are not similar, and generate alike networks. Therefore, we consider that $0.0004\lesssim p\lesssim 0.1$ is an interesting range of edge densities because only the most similar pairs of time series are connected and the networks are not fragmented into too many components. This is also well in line with the settings used in the vast majority of previous climate network studies. For a more in-depth analysis, we have selected four values of $p$ based on its influence on the network topology: $0.001,0.01,0.05$, and $0.1$. The results of our following analyses will be reported for those four cases.

3.2 Global network features

In the following, we further analyze a few selected network properties that are commonly used to classify complex network into certain categories like scale-free and small-world networks [16]. Figure 4A shows the values of the smallworldness index [59] for the different distance functions and four values of the edge density. According to its definition, values of the smallworldness index larger than one (dashed line) are indicative of the small-world feature (i.e., high $C$ along with low $L$). The exact values for the network measures can be found in Appendix B, Tab. 5, 6, 7 and 8. Our results show that all generated networks indeed exhibit the small-world property, except for the correlation, cross-correlation, mutual information, MIC, STS and CORT functions when $p=0.001$. These results are in agreement with previous studies reporting a small-world effect in climate networks [4, 7]. Our results show that this effect is also present when other distance functions are considered. On the one hand, the high clustering coefficient results from the spatial continuity of the underlying physical fields [4] which tends to connect close cells (nodes) and thereby form triangles. On the other hand, the existence of long-range connections (teleconnections) explains the small average path length. These long distance edges are responsible for an efficient propagation of local perturbations through the whole network, as will be further discussed in Section 3.4. We emphasize again that one reason for the small-world effect to emerge across all climate networks is the spatial embedding of the nodes. Although we found evidences of small-worldness for almost all distance functions, some or all of them might be biased, as demonstrated by [61]. Studying whether or not this already explains the observed characteristics completely would be possible by using sophisticated spatial network surrogates [23], the consideration of which is however beyond the scope of the present work.

Small-world networks like those presented by the studied climate networks can be further classified according to their degree distributions into three types: scale-free, broad-scale and single-scale [62]. Scale-free networks have a power law distribution, broad-scale networks have a power law regime followed by a sharp cutoff, and single-scale networks present a faster than algebraic decay, e.g., exponential or Gaussian. In Fig. 4B and C, we show the cumulative degree distributions for all networks separated by edge densities in a log-log and linear-log scales, respectively. Our results indicate that all distance functions generate broad-scale networks for small edges densities ($p=0.001$), i.e., they present a scale-free interval followed by sharp decays. For larger edge densities, all distance functions generate single-scale networks.

3.3 Network similarities

We next investigate the similarities in the network structures created using different time series distance functions. As in our previous analyses, we consider four edge densities: $0.001,0.01,0.05$, and $0.1$. For each value of $p$, we construct the networks using all 29 distance functions (Tab 1) and compare the resulting networks using the Hamming distance [58] (Eq. 3). Therefrom, we build a distance matrix $\mathbf{D}^{H}$ whose element $D^{H}_{A,B}$ is the Hamming distance $d_{H}(A,B)$ between the networks $A$ and $B$. Then, we use hierarchical agglomerative clustering (using Ward’s method) to cluster the different networks.

The obtained clustering results are shown in Fig. 5. The distance matrices for the four edge densities can be found in Appendix C (Fig. A). Our results suggest that, for all four edge densities, all PDF-based distance functions (Tab. 1) except the mutual information form one cluster (blue dotted lines in Fig. 5). This indicates that all of them generate similar networks independent of $p$. In the literature, a large number of studies has constructed functional climate networks using Pearson correlation with or without lag, while a minority has used mutual information (MI) [4, 5, 7, 9, 10]. According to our results, these three approaches generate similar networks (red dotted lines in Fig. 5), which confirms previous results [4, 63]. We also find that the Pearson correlation and cross-correlation generate similar networks, which suggests that previous studies using both measures are compatible and that the lag does not markedly affect the network topology.

Figure 5 also presents a comparison between the dendrograms and their cophenetic correlations. The idea here is to verify if the clustering hierarchies are consistent for different edge densities. The black lines between dendrograms indicate the positions of one network in different dendrograms. Strong entanglement indicates large differences between the dendrograms. Our results show that there exists some entanglement between $p=0.001$, $0.01$ and $0.05$ (Fig. 5A and B) that is mainly caused by the small number of links in the networks. Even though, the clustering hierarchies can be considered similar since their cophenetic correlations are high ($>0.89$), especially at higher edge densities (Fig. 5C).

3.4 Teleconnections

We have already shown that networks constructed with different distance functions exhibit distinct topologies. We have also observed that almost all resulting networks display the small-world feature. If the networks are small-world and have different topology, the long-distance edges responsible for the small-world feature should be different. This indicates that distinct distance functions can capture different long-distance relationships which can be caused by different underlying climate effect.

In the following, we further analyze the differences between the topologies generated by different distance functions. As described in Sec. 2.3, we define teleconnections as the most similar pairs of time series whose grid cells have great-circle distances longer than $5000$ km [60]. Atmospheric teleconnection patterns can be generated by both, internal atmospheric processes and forcing from surface conditions. According to [64], such teleconnections enable the atmosphere to act like a "bridge" between different parts of the ocean and while enabling the ocean to act like a "tunnel" linking different atmospheric regions. These long-distance links are of great importance for climate science because they can provide significant skill to statistical weather forecasts from synoptic to seasonal scales and even beyond [27, 28].

Our methodology (Sec. 2.3) consisted on clustering the network topologies and extracting the centroid network of each of cluster. Figure LABEL:fig:teleconnections shows the 500 strongest teleconnection related links for the eight different centroid networks. As already demonstrated in previous works [4, 10, 63], the correlation captures especially the influence of the El Niño-Southern Oscillation (ENSO). ENSO is a very important teleconnection pattern responsible for climate variability all over the globe [24, 25, 26]. This phenomenon originates from nonlinear ocean-atmosphere interactions in the tropical Pacific ocean connecting large-scale oceanic sea surface temperature (SST) anomalies (El Niño) with the large-scale atmospheric Southern Oscillation (SO), which is characterized by a strong sea-level pressure seesaw between French Polynesia and North Australia [65, 66]. This pattern is also captured by Euclidean, DTW, INTPER, and ACF distances. The teleconnections found by CDM and Gower distances do not exhibit coherent spatial patterns, which is expected in physical mechanisms [28]. This lack of coherence does not imply that patterns do not exist but that we cannot identify any using our methodology.

The Euclidean, DTW, and PER distances also capture teleconnections similar to those associated with the Pacific-South America (PSA) pattern and the Antarctic Oscillation (AAO). The PSA modes are unique features of atmospheric variability in the Southern Hemisphere in the subpolar/polar regions of the South Pacific. The PSA pattern was suggested as part of a stationary Rossby wave train, which is usually generated by changes in tropical convection [67]. According to [68], the PSA pattern is related to sea surface temperature (SST) variability over the central and eastern Pacific at decadal scales, and it is a response to the El Niño–-Southern Oscillation (ENSO) at interannual time scales. The PSA can also be associated with the quasi-biennial component of ENSO, and the strongest connections occur during austral spring. This pattern represents a zonal wave with wavenumber 3 and a well-defined wave train from the tropical Pacific and the Indian Ocean to South America with a large amplitude in the PSA sector [68]. The zonal wavenumber represents the number of waves at a given latitude. The large amplitude in the PSA sector means that the cold air advance over a region, with impacts in the regional climate [69]. PSA is also related to the AAO, which is characterized by zonally symmetric structures with opposite signs between the Antarctic region and latitudes of about 45$\degree$S [70]. The annular pattern of AAO can be considered as a presence of anomalies with latitudinal variation between the regions of extratropical and subtropical latitudes of the Southern Hemisphere [71]. The positive phase occurs when there is a predominance of negative SST anomalies close to Antarctica, and positive in the mid-latitudes, resulting in a shift of extratropical cyclones and cold fronts towards Antarctica [72].

Through INTPER, it is possible to observe the relationship between the eastern Pacific and tropical Atlantic temperatures that can be explained by the Southern Oscillation [73]. The ACF and INTPER show an interesting inter-basin connection between the Indian, Pacific and Atlantic oceans. This connection can be explained by the Indian Ocean Dipole (IOD). The positive phase of the IOD is associated with warmer (colder) than normal SST off the coast of East Africa (Java-Sumatra) [74]. The signal of the associated atmospheric teleconnection can be noticed far from the source region [75, 76, 77], also affecting the SSTs in the Atlantic Ocean near South America [74]. Empirical orthogonal function (EOF) analysis indicates that 30% of the total variation of anomalous Indian Ocean SSTs is during ENSO events and the dipole mode (IOD) explains about 12% of the associated variability [78].

Our results further show that different distance functions are capable to capture different patterns and relationships that are not very well understood, such as those shown in ACF and INTPER based networks. These patterns represent an enhancement of existing knowledge and should be further investigated. Furthermore, some distance functions also indicate interesting interactions in the meridional direction, connecting the two hemispheres. Those interactions are not trivial as the atmosphere and ocean tend to flow zonally as constrained by the conservation of angular momentum [64], which suggest that studies to explore these connections can be made through the presented functions.

4 Conclusions

In this work, we have analyzed the effect of 29 different time series distance functions on the topology of functional climate networks. We have started by studying how the edge density $p$ affects the resulting network topology. When $p\lesssim 0.0004$, only a few very similar nodes are connected and the networks present many components. On the other hand, when $p\gtrsim 0.01$ the networks become densely connected and even not very similar nodes are mutually linked. The interval $0.0004\lesssim p\lesssim 0.1$ manifests a transition phase where the networks get more and more connected. We consider this interval as relevant for the study of climate networks, at least for the specific dataset employed in our work. The resulting networks exhibit a (likely geometrically induced) small-world effect and have a degree distribution that decays faster than a power law distribution.

A large body of previous studies have used the Pearson correlation coefficient, covariance or cross-correlation to measure similarity between time series. We have shown here that these measures generate similar networks. Our results also demonstrate that other distance functions are capable of generating different network topologies. As all these networks exhibit similar small-world characteristics, yet with distinct topologies, the long-range connections responsible for this effect tend to be different.

It is important to remember that, although we found evidences of small-worldness for almost all distance functions, this effect may be biased by the network construction process [61], which represents a limitation in the methodology here proposed. Nevertheless, we have observed that distinct distance functions can capture different teleconnection patterns. These teleconnections play an import role in the climate system since they generally represent energy transfer between remote regions [26] and can be used for improving statistical weather forecasts at various time scales [27, 28]. Although the structure of several teleconnections patterns has been known, the reasons why they emerge are not always well understood [25].

One advantage of functional climate networks is the possibility to inspect the climate dynamics from a complementary perspective. We have shown here that alternative distance functions capture different teleconnections that present patterns similar to some known climate phenomena like the Pacific-South America (PSA) pattern, Antarctic Oscillation (AAO), and Indian Ocean Dipole (IOD). Phenomena like these have not been investigated in previous works using climate networks. A few distance functions, such as CDM and Gower (Fig. LABEL:fig:teleconnections), did not present coherent spatial patterns. It does not mean that a pattern does not exist, but that we cannot identify any using our methodology. We believe that this is not a limitation in our methodology, but an indicative that some distance functions can be inappropriate for the study of climate. To conclude, we believe that our results provide a new perspective for climate science because they allow the investigation of climate patterns from a different point of view. This advanced approach could finally lead to a better understanding of how global climate patterns emerge.

In summary, we have observed differences in the topologies of climate networks generated by different distance functions. These distance functions capture non-trivial climatological phenomena that previous studies on climate network have not found. This research opens many opportunities for future works. One obvious next step is the thorough evaluation of the climate patterns and the deep reasons behind the topological differences. This analysis can be easily extended to different climate variables using the distance functions used in this study. Other distance functions may also be considered, for example, synchronization-based measures and causality analysis can be applied. Ultimately, taking temporal variations in the corresponding climate patterns into account, such an investigation can also be extended into a temporal network framework.

Acknowledgments

This research received financial support from the São Paulo Research Foundation (FAPESP) grants 2017/05831-9 and 2015/50122-0. The authors also acknowledge the National Council for Scientific and Technological Development (CNPq) for its financial support. This research was developed using computational resources from the Center for Mathematical Sciences Applied to Industry (CeMEAI) funded by FAPESP (grant 2013/07375-0). RVD has received funding by the German Federal Ministry for Education and Research (BMBF) via the BMBF Young Investigators Group CoSy-CC2 (Complex Systems Approaches to Understanding Causes and Consequences of Past, Present and Future Climate Change, grant no. 01LN1306A), the Belmont Forum / JPI Climate project GOTHAM (Globally Observed Teleconnections and Their Representation in Hierarchies of Atmospheric Models, grant no. 01LP16MA) and the JPI Climate / JPI Oceans project ROADMAP (The Role of Ocean Dynamics and Ocean-Atmosphere Interactions in Driving ClimAte Variations and Future Projections of Impact-Relevant Extreme Events, grant no. 01LP2002B).

Computer Code Availability

The method and the experiments in this paper were implemented using the R programming language [79] with the packages igraph [80], TSclust [81], and TSdist [82]. The R software, the packages, and our implementation are all open-source (GPL) and freely available for download. Our source codes can be obtained at https://lnferreira.github.io/climate_networks_R.

Appendix A Definitions of time series distance functions