Global teleconnectivity structures of the El Ni\~no-Southern Oscillation and large volcanic eruptions -- An evolving network perspective

Recent work has provided ample evidence that global climate dynamics at time-scales between multiple weeks and several years can be severely affected by the episodic occurrence of both, internal (climatic) and external (non-climatic) perturbations. Here, we aim to improve our understanding on how regional to local disruptions of the"normal"state of the global surface air temperature field affect the corresponding global teleconnectivity structure. Specifically, we present an approach to quantify teleconnectivity based on different characteristics of functional climate network analysis. Subsequently, we apply this framework to study the impacts of different phases of the El Ni\~no-Southern Oscillation (ENSO) as well as the three largest volcanic eruptions since the mid 20th century on the dominating spatiotemporal co-variability patterns of daily surface air temperatures. Our results confirm the existence of global effects of ENSO which result in episodic breakdowns of the hierarchical organization of the global temperature field. This is associated with the emergence of strong teleconnections. At more regional scales, similar effects are found after major volcanic eruptions. Taken together, the resulting time-dependent patterns of network connectivity allow a tracing of the spatial extents of the dominating effects of both types of climate disruptions. We discuss possible links between these observations and general aspects of atmospheric circulation.


Introduction
The empirical analysis of climate data is fundamental to understand the evolution and develop more accurate methods for forecasting of climate phenomena like El Niño. Typically, such data sets comprise time series representing temperature, precipitation or other climate variables observed at distinct locations distributed around the globe. Their common properties include long-range spatial and often also temporal correlations (Fraedrich and Blender, 2003), interactions at and among multiple scales (Paluš, 2014) and nonlinearity (Dijkstra, 2013). With the Earth's surface being subdivided into regions for which individual "grid points" and associated localized records of climate variability are considered representative, the evolution of the climate system can be approximately described by a high-dimensional multivariate time series composed of a multitude of interdependent signals.
While the analysis of such big climate data sets has been traditionally attempted by means of classical statistical approaches like empirical orthogonal function or maximum covariance analysis (von Storch and Zwiers, 2003), it has recently been realized that these methods exhibit fundamental intrinsic limitations, including their linearity and associated condition of pairwise orthogonal patterns (Gámez et al., 2004). As a consequence, the traditional view that the corresponding decompositions of global spatio-temporal covariability patterns actually provide dynamical (or at least statistical) modes that unambiguously coincide with specific key climatic processes has been abandoned (Monahan et al., 2009). Taken together, there is growing evidence that the application of traditional linear methods of signal processing and the climatic interpretation of their results are severely affected by the dynamical complexity of the involved processes.
During the last years, complex network representations of climate variability have been developed (Tsonis et al., 2006;Donges et al., 2009a, b;Tsonis et al., 2011;Tsonis and Swanson, 2012;Steinhaeuser et al., 2012;Peron et al., 2014;Ciemer et al., 2017) and demonstrated to provide a suitable approach for relieving some of the aforementioned concerns (Donges et al., 2015b). In this nonlinear statistical framework, referred to as functional climate network analysis, the individual grid points or cells are considered as nodes of a spatially embedded graph. Connections among these nodes are established according to similarities between the individual (local) climate time series (Tsonis et al., 2006;Donges et al., 2009b;Donner et al., 2017). By construction, the network structures thus obtained highlight essential statistical interrelationships among spatio-temporal climate data (Donges et al., 2009b).
The application of functional climate networks has already provided several important insights. For instance, centrality measures, such as betweenness centrality, have been found to serve as tracers of global circulation patterns in the atmosphere and oceans (Donges et al., 2009a). Moreover, climate networks have been used to identify dipole patterns which represent pressure anomalies of opposite polarity appearing in two different regions simultaneously (Kawale et al., 2011). The study of the coupling structure between interdependent climate variables (Donges et al., 2011), the temporal evolution and teleconnections of the North Atlantic Oscillation (Guez et al., 2012(Guez et al., , 2013, the distinction of different types of El Niño phases (Radebach et al., 2013;Wiedermann et al., 2016) and the prediction of the latter (Ludescher et al., 2013(Ludescher et al., , 2014 have also been subjects of corresponding recent investigations. Many of the aforementioned methodological achievements have been integrated in open source software packages (Donges et al., 2015a) contributing to the increasing use of functional network analysis in climatological studies (Donner et al., 2017).
One rather fundamental property of large networks is their (possibly hierarchical) organization in terms of communities -an aspect that has also been addressed recently in the context of key patterns in climate data (Tsonis et al., 2011). Here, a community is a subset of densely connected nodes which exhibit only few interactions with the rest of the network (Newman, 2006b;Fortunato, 2010). In a climate network context, communities would ideally have some climatological interpretation. Specifically, Tsonis et al. (2011) argued that each community in a climate network should be considered as a subsystem which operates relatively independent of the other communities. Besides corresponding connectivity structures in individual climate variables, community detection algorithms (Fortunato, 2010) can also be used to detect multi-variable clusters (Steinhaeuser et al., 2010).
In this paper, we analyze global surface air temperature data in terms of functional climate networks and demonstrate the close relationship between El Niño and La Niño episodes as well as strong volcanic eruptions on the one hand, and temporal changes in the modular organization of the resulting networks on the other hand. For this purpose, we study 3.4

ENSO-big
Agung El Chichon Pinatubo Figure 1. Main regions of interest used within this paper. Sets of blue dots labelled with "El Chichon", "Agung" and "Pinatubo" indicate grid points within a 5 • radius around the corresponding volcanoes. The numbers 3, 3.4 and 4 identify the corresponding Nino regions (cf. Tab. 1) commonly used for defining characteristic indices of ENSO variability. The region "ENSO-big" will be removed from the complete global data set when analyzing the spatial imprints of volcanic eruptions to ensure that ENSO-related effects are excluded.
the teleconnectivity structure in the climate system in terms of spatial fields of two network properties that represent the number of strong statistical connections, as well as the average spatial distance between the connected grid points. In addition, the associated temporal variations are traced by some scalar-valued global network characteristics.
The remainder of this paper is organized as follows: Sect. 2 provides brief information on the climatological background of ENSO and volcanic eruptions as the two types of major climatic disruptions studied in this work. The data and methodology used in this work are described in detail in Sect. 3. Finally, our results are presented and discussed in Sect. 4, followed by concluding remarks.

El Niño-Southern Oscillation
Among the dominant teleconnectivity patterns in the global climate system, the El Niño-Southern Oscillation (ENSO) is the probably most remarkable phenomenon in terms of both, the magnitude of associated variations in sea-surface temperature (SST) and sea-level pressure, as well as its specific impacts on different aspects of regional climate variability worldwide (Trenberth, 1997). During the positive phase (El Niño) of this complex oscillation of the coupled atmosphereocean system in the tropical Pacific, the eastern tropical Pacific exhibits some anomalous warming with respect to "normal" mean conditions, while the negative phase (La Niña) is characterized by a corresponding cooling. In comparison with the normal mean climatology, this surface temperature region latitudes longitudes Nino1+2 Table 1. Overview on different regions commonly used for defining characteristic temperature-based indices associated with ENSO variability. In addition, we include the definition of the "ENSO-big" region studied in this work, which corresponds to the region that is discarded in our analyses of the impacts of strong volcanic eruptions on global temperature teleconnectivity.
anomaly results in marked shifts of key atmospheric pressure systems, modifying the large-scale circulation and, thus, leading to prominent shifts of, e.g., precipitation patterns. It has been shown that effects of both ENSO phases can be observed in remote regions including North and South America, Africa, the Indian subcontinent, and even Antarctica (Ropelewski and Halpert, 1987;Dai and Wigley, 2000;Neelin et al., 2003;Turner, 2004;Clarke, 2008;Sarachik and Cane, 2010).
The long-term variability of ENSO is characterized by some irregular oscillations with a period of 2 to 7 years and remarkable variations in the associated characteristic frequencies and amplitudes of both, temperature and pressure anomalies. Following its prominent spatial structures in tropical SST and sea level pressure, ENSO is commonly traced by indices that take up the variability of the aforementioned observables in some key region of the tropical Pacific ocean. Notably, a set of indices has been defined in terms of average SST anomalies taken over distinct regions in the eastern and central tropical Pacific, referred to as Nino1+2, Nino3, Nino4 and Nino3.4, respectively (Trenberth and Stepaniak, 2001) (see Fig. 1 and Tab. 1). In this work, we will utilize the socalled Ocean Niño Index (ONI) for differentiating between different phases of ENSO. It is defined as the running threemonth mean SST anomaly for the Niño 3.4 region (5 • N-5 • S, 120 • -170 • W) in comparison with centered 30-year base periods that are updated every 5 years (NCEP, 2017). When the ONI exceeds 0.5 • C for at least five consecutive months, the corresponding situation is classified as an El Niño, and the magnitude of the ONI is taken as an indicator of the strength of the corresponding event. In turn, if the ONI drops below −0.5 • C for at least five consecutive months, this indicates a La Niña episode.
In the last years, it has been recognized that the commonly observed spatial patterns associated with El Niño (as well as La Niña) related SST anomalies are far from being homogeneous across the set of observed events. Consequently, it has been suggested to further distinguish both phases into two respective flavours (see Wiedermann et al., 2016, and references therein). The first type is the classic or East Pa-
cific (EP) El Niño (Rasmusson and Carpenter, 1982;Harrison and Larkin, 1998), which is localized in the eastern tropical Pacific and characterized by strong positive SST anomalies close to the western coast of South America. Opposed to this, the El Niño Modoki or Central Pacific (CP) El Niño exhibits marked SST anomalies in the central tropical Pacific close to the dateline. Notably, both spatial structures (EP and CP) can be observed in the context of La Niña, too. Noticing that there have been contradictory classifications in the literature for some past ENSO events, Wiedermann et al. (2016) recently presented a new indicator for the ENSO flavor based on functional climate networks. In the remainder of this paper, we will follow their classification, which is summarized in Tab. 2.

Volcanic eruptions
Besides distinct ENSO episodes and their known global climate impacts, another type of events that can substantially affect climate at large spatial and temporal scales are strong volcanic eruptions. Similar to El Niño and La Niña episodes, such events can result in large-scale spatially coherent cooling trends due to modifications of the radiation balance by changes in atmospheric chemistry and the shielding effect of volcanic aerosols in the stratosphere. Subsequently, such cooling can again cause changes of precipitation and temperature patterns from synoptic (weather) time scales to relatively persistent multi-annual effects (Robock, 2000) and even trigger long-lasting climate disruptions (Miller et al., 2012). In the past decades, several large volcanic eruptions have injected up to some millions of tons of sulfur dioxide into the atmosphere, which can get rapidly distributed around the globe once entered the stratosphere.
In this study, we focus on the global effects of the three major volcanic eruptions during the second half of the 20th century. Within this period, the largest and most influential event, the Mount Pinatubo eruption (McCormick et al., 1995), took place between April and September 1991 in the Philippines, followed by the Mount Agung eruption in Indonesia (February 1963to January 1964 and the El Chichon eruption (March to September 1982) in Mexico (see Fig. 1).

Data
In this study, we use daily mean surface air temperature (SAT) data (at sigma level σ = 0.995) from the National Center for Environmental Prediction (NCEP) and National Center for Atmospheric Research (NCAR) Reanalysis I project (Kalnay et al., 1996;Kistler et al., 2001): The data cover the years 1948-2015 at a global grid with equi-angular spatial resolution of 2.5 • in both latitude and longitude, thus comprising 10, 512 individual temperature time series. Note that we found the inclusion of the land areas important, hence we chose SAT instead of the aforementioned SST. In order to remove leading order effects of seasonality in the temperature recordings, the long-term average temperatures for each calendar day of the year have been subtracted from the raw data independently for each grid point, resulting in so-called SAT anomalies.
Equi-angular gridded data have, by construction, a higher density of grid points at the poles than around the equator, which would result in systematic biases of statistical characteristics overemphasizing the polar regions with apparently more data if not properly accounted for. For the latter purpose, area-weighted measures have been developed and subsequently applied in recent works (Tsonis et al., 2006;Heitzig et al., 2012;Wiedermann et al., 2013). As an alternative, we follow here the approach of Radebach et al. (2013), where the original data have been remapped onto a grid with a much higher spatial homogeneity. Specifically, we use an icosahedral grid as described by Heikes and Randall (1994), which finally leads to a decomposition of the Earth's surface into Voronoi cells of almost the same area. In the present case, the corresponding remapping procedure results in a set of N =10,242 grid points that exhibit a narrowly peaked distribution of geodesic distances between direct neighbors.
In Radebach et al. (2013), the time series associated with each new grid point have been determined based upon the values at the respective four surrounding grid points in the original equi-angular grid. In this work, we use a slightly different approach by taking the four closest points in space instead, which in some cases may deviate from the former setting. This modification is motivated by the fact that the consideration of the spatially closest "observational" values may provide a better approximation of climate variability at the new grid point. Moreover, these spatial neighbours can be determined efficiently using spatial search trees. Due to the commonly rather large spatial correlation length of the SAT field (as compared to other climate variables like precipita-tion) and its resulting spatial smoothness, we do not expect the time series resulting from both algorithmic variants to differ markedly.
Finally, we note that when using the global data set as described above, the temporal correlations associated with the key ENSO region and the surrounding parts of the Pacific ocean are known to dominate climate variability globally. This leads to undesired outcomes when aiming to resolve the effects of individual volcanic eruptions on global temperature patterns, since they might be masked by ENSO variability, especially in cases where the corresponding effects take place simultaneously with some El Niño (or La Niña) event. In fact, strong tropical volcanic eruptions have even been suggested to serve as triggers for El Niño phases (Khodri et al., 2017).
In order to account for the problem of temporal cooccurrence between the effects of volcanic eruptions and ENSO events, we are going to use the full set of data when studying the effects of ENSO on global temperature teleconnectivity, while excluding the main ENSO region and its surroundings (referred to as "ENSO-big" in Fig. 1 and Tab. 1) when studying the impacts of specific volcanic eruptions. Note that this excluded region has been chosen rather large on purpose (as an outcome of more systematic studies with variable regions to be discarded, which are not further discussed here for brevity) such as to ensure an as complete as possible separation between the direct ENSO impacts and the effects of volcanic eruptions, especially in cases of simultaneous events. In fact, when considering the full global SAT data set in the context of the impacts of volcanic eruptions, only the signatures of the Mount Pinatubo eruption are clearly visible (Radebach et al., 2013). Alternative strategies for reducing the impact of ENSO variability in order to highlight the climate effects of other phenomena like volcanic eruptions might include conditioning out the effect of ONI or other representatives of the ENSO state on the local SAT variability at each grid point prior to network analysis. We outline further investigations that make use of such approaches as a subject of future research.

Functional climate network analysis
Functional climate networks provide a coarse-grained spatial representation of the co-variability structure among globally or regionally distributed measurements of some climate variables (Tsonis et al., 2006;Tsonis and Swanson, 2012;Donner et al., 2017). Starting from a set of records of the variable of interest, the geographical positions associated with the individual time series are identified with the N nodes of some abstract network embedded on the Earth's surface. The connectivity of this network is then formed by establishing links between pairs of these nodes that fulfill some statistical similarity criterion (see below). Thus, links in such climate networks represent strong statistical associations between climate variability at different points in space. These associa-tions may potentially indicate certain climatic processes or mechanisms interlinking the variability at the corresponding locations. Hence, the resulting linkage structure is referred to as functional connectivity.
Like other undirected and unweighted networks, functional climate networks are conveniently represented in terms of their adjacency matrix A = (A ij ), where A ij = 1 indicates the existence of a link between node i and node j, while A ij = 0 corresponds to an absence of such a link. In our specific case, the matrix A is time-dependent, since the spatial co-variability structure of the SAT field changes with time. In such a case, we speak of an evolving climate network (Radebach et al., 2013).

Network generation
According to our considerations presented above, we take the grid points of the icosahedral grid constructed by remapping the original NCEP/NCAR reanalysis data as nodes of an evolving SAT network (i.e., we consider a fixed node set that does not change over time).
For establishing the time-dependent link set, we consider sliding windows in time covering a set of days {d} = [d 0 , d 0 + ∆d] with a width of ∆d = 365 days and mutual offset of 183 days between subsequent windows. Each of these windows is labeled with the corresponding midpoint date d mid = d 0 + ∆d/2. Our choice of 1-year windows is motivated by the fact that seasonality may not be not completely accounted for when using shorter windows, since the consideration of SAT anomalies as defined above does not exclude seasonality in the local higher-order statistical characteristics of SAT after correcting only for the mean climatology. Moreover, due to different distributions of land and water masses on both hemispheres of the Earth (with difference persistence properties of SAT), the resulting spatial covariability structure manifested in the climate network topology may undergo seasonal variations as well, which could affect the results of our analysis presented in this work.
From the seasonally adjusted temperature data at each grid point during a time window (corrected for the window-wise mean), T i ({d}) (with i denoting the respective grid point), we compute the matrix of pairwise Pearson's correlation coefficients where • {d} and σ {d} (•) denote the mean value and standard deviation of the respective variable taken over the time window {d}. From this matrix, we identify the entries (i.e., pairs of nodes) with the highest absolute values of mutual correlations |c ij |. Specifically, in this work, we consider the 0.5% strongest pairwise statistical similarities among all nodes per window, i.e., where Θ(•) is the Heaviside function, q |c|,0.995 (d mid ) is the 99.5-percentile of the distribution of absolute correlation values for the time window centered at d mid , and δ ij denotes Kronecker's Delta. A (d mid ) is now the mathematical representation of our evolving climate network. In the following, we omit the explicit time dependence of A (also in all network properties) for brevity. Note that according to our construction, A is symmetric.

Node degree
The degree k i of a node i is defined as the number of links connected to i, It represents how densely a node is connected within the network. In case of a functional climate network, the degree can thus be considered as a proxy for the importance (or centrality) of a certain grid point in the spatio-temporal correlation structure of the variable of interest.
In the following, we refer to network measures like the degree, which provide a characteristic value specific to each individual node i, as local network characteristics. We call the full set of their values taken together with the associated spatial positions of all nodes a field.

Average link distance
Another local network characteristic that defines a field of values upon a functional climate network is the average link distance (Donner et al., 2017) of a node i, which is defined as with d ij being the normalized spatial distance between two nodes i and j. As a proper normalization, we choose here the largest possible (shortest) distance between two points on the Earth's surface, i.e., half of the circumference of the Earth, u Earth , so that d ij = 2D ij /u Earth where D ij is the geodesic distance between nodes i and j. A low average link distance indicates that i has very localized connections, while a high value points to a node with long-ranging spatial connections. This measure is closely related with the total distance of a node with respect to the rest of the network as previously used by  in the context of functional climate networks. We emphasize that the average link distance must not be confused with the conceptually related average path length, where d ij , as defined above, would be replaced by the minimum number l ij of links separating two nodes i and j in the network (i.e., where i and j do not need to be directly connected).
Taking the average of d i over all nodes i of the network gives the global average link distance

Transitivity
The network transitivity quantifies how strongly the connectivity of a network is clustered. Mathematically, it describes the degree to which the network's adjacency property is transitive, i.e., the fraction of cases in which the presence of two links between nodes i and j as well as i and k is accompanied by a third link between j and k. Mathematically, this is expressed as (Boccaletti et al., 2006;Radebach et al., 2013) Like the global average link distance d (but unlike degree and average link distance), T does not define a field, but returns one single single scalar value for each network.

Modularity
The concept of modularity was introduced into network science by Newman and Girvan (2004) to measure the degree of heterogeneity within the network structure, i.e., how well different groups of nodes can be distinguished that are densely connected within each group, but only sparsely among different groups. In the case of a climate network, modularity provides a single scalar-valued characteristic property that discriminates between a relatively homogeneous link placement (low modularity) and the existence of (commonly regionally confined) clusters of nodes (time series) that exhibit relatively coherent variability (high modularity). The definition of modularity relies upon a partitioning of the network into meaningful subgraphs. Up to a multiplicative constant, it counts how many links are clustered within these subgraphs and compares this value with the expected number of links inside these subgraphs if the network were random, where m is the total number of links and ∆ ij an indicator function informing whether or not two nodes i and j belong to the same subgraph in the considered partition. The individual subgraphs in the partition that maximizes the modularity Q are called communities. The higher the modularity, the more split-up (or modular) the network. Accordingly, community detection by modularity maximization has become a common tool for cluster analysis.
While the above definition of modularity is mathematically precise, its maximization is a hard computational problem and can only be achieved by making use of suitable heuristics. Various estimation algorithms have been proposed (Fortunato, 2010). It should be emphasized that many of them can result in suboptimal solutions. Thus, a good choice of the algorithm is important for obtaining reliable results. In this work, we employ the WalkTrap method introduced by Pons and Latapy (2006). By comparing the results provided by this algorithm with those of other methods (cf. Appendix A), we have found that the WalkTrap solution exhibits comparably high values of modularity and relatively stable values in case of strongly overlapping windows (i.e., in cases where the considered data do not change much).

Regionalization of field measures
As detailed above, node degree and average link distance constitute two important local network characteristics. In some of our following investigations, it will be useful to study the associated spatial fields in full detail. However, when focusing on the specific impacts of certain climate phenomena, it can be beneficial to perform a regionalization of these measures. Specifically, for a subset of nodes X ⊆ {1, . . . , N } representing a certain part of the globe, a regionalized version of the degree would be given as where |X | denotes the number of nodes in the considered set. As a consequence, we can not only assign a degree value to an individual node, but also (as a mean degree) to a subgraph. Note that this regionalized degree differs from the concepts of cross-degree and cross-link density between subgraphs (Donges et al., 2011), since unlike k X , the latter exclude contributions due to links between nodes within X in their definition. For the average link distance, the corresponding regionalized property d X can be defined in full analogy.
Below, we detail some reasonable choices for X to be utilized in the context of the present work, which focus on specific spatially contiguous regions of the Earth's surface that are associated with ENSO or volcanoes with strong past eruptions.

El Niño-Southern Oscillation regions
As already detailed in Sect. 2.1, there exist a variety of indices that measure the "strength" of a particular ENSO state. Among others, four Nino regions (1+2, 3, 4 and 3.4) have been defined to capture SST anomalies in different parts of the tropical Pacific.
The regionalization approach described above can be applied to these four regions by taking all nodes located within the respective spatial domains and apply Eq. (8). Thereby, we obtain a set of eight new scalar-valued characteristics: k Nino1+2 , d Nino1+2 , k Nino3 , d Nino3 , k Nino4 , d Nino4 , k Nino3.4 and d Nino3.4 . In order to reduce this vast amount of information, in what follows, we will not make use of the (anyway less frequently studied) Nino1+2 region, but focus on the Nino3.4 region (which is also the basis of the nowadays most common ONI index) and its two contributors, Nino3 and Nino4.

Volcano regions
The locations of the three volcanoes responsible for the largest eruptions of the recent decades are shown in Fig. 1. To obtain interpretable information on the (tele-) connectivity induced by these eruptions, we need to integrate the connectivity properties of a sufficiently large amount of meaningfully chosen grid cells. As a first attempt, we therefore take the area within a radius of 5 • around each volcano as basis for the regionalization procedure of k i . This leads to the three observables k Pinatubo , k Agung and k Chichon . For the average link distance, one could again proceed in a similar way.
However, the aforementioned choice might not be optimal, since symmetric spatial regions in the near-field do not necessarily exhibit the strongest persistent temperature effects after an eruption. Instead, the specific local meteorological conditions (especially wind fields) during the eruption period largely control the three-dimensional patterns of atmospheric aerosol concentrations and, hence, the position of the strongest mid-term cooling to be expected. Accordingly, the induced teleconnectivity can be more evident within regions that have been shifted with respect to the locations of the volcanoes. To account for this, we also calculate regionalized degrees for accordingly shifted regions (see Sect. 4.2 for details), denoted as k Pinatubo , k Agung and k Chichon . Here, the specific regions have been selected according to an examination of the resulting degree fields of the SAT networks for time windows succeeding the individual eruptions and the corresponding wind fields, seeking for the timing and position of the strongest anomalies in the degree field that could be attributed to each eruption (see below). Note that although a volcanic eruption may start relatively abruptly, its larger-scale atmospheric effects commonly become effective only with a considerable delay of several months or more (Robock, 2000;McCormick et al., 1995).

Results
In the following, we present the results of our functional network analysis of global SAT patterns with a focus on the associated imprints of ENSO. Subsequently, we turn to analyzing and discussing the excess connectivity induced by strong volcanic eruptions.

El Niño-Southern Oscillation
Let us start with investigating the global effects of ENSO on the spatio-temporal co-variability structure of global SAT. From a complex network perspective, this problem has already been addressed in a variety of previous studies (e.g. Yamasaki et al., 2008;Gozolchiani et al., 2008;Radebach et al., 2013;Wiedermann et al., 2016;Fan et al., 2017, and various others), making use of different approaches for constructing network structures from global climate data. However, none of these works has considered the complementarity between topological and spatial network properties in great detail, nor utilized the concepts of modularity and global average link distance that constitute key aspects of this paper and provide important new insights as demonstrated in the following.

Global network properties
The network transitivity T has been shown by Wiedermann et al. (2016) to systematically discriminate between the EP and CP flavours of both, El Niño and La Niña. While this reference used an area-weighted version of T and included information on the total pairwise correlation strength instead of just binary adjacency information, we follow here the approach of Radebach et al. (2013) in using remapping onto an icosahedral grid. Figure 2a shows the corresponding results obtained using our slightly modified data set, which are qualitatively almost indistinguishable from those of the two aforementioned studies as expected 1 .
In order to further quantify the strength of teleconnectivity in the global SAT field, the network modularity Q provides a prospective candidate measure that has not yet been exploited for this purpose in previous studies. Recall that a high modularity indicates a fragmented network, whereas low values would point to a relatively homogeneous connectivity structure of the network as a whole. Hence, a marked decrease in modularity could indicate an increase in the degree of organization of the global SAT network, i.e., a tendency towards more balanced co-variability in global temperatures. Figure 2b shows that most time intervals that are characterized by elevated values of network transitivity actually exhibit a marked reduction in modularity. Consistent with previous findings of Radebach et al. (2013), most of these time windows in fact coincide with either some El Niño or La Niña phase, indicating again the global impact of these episodes in terms of equilibrating spatial co-variability in the Earth's SAT field. This can be considered as an expected signature of emerging teleconnectivity. Note that taken alone, this process would not necessarily imply a stronger synchronization (as studied by, e.g., Maraun and Kurths (2005)) be-1 Note, however, that Fig. 2 1 9 5 0 1 9 5 5 1 9 6 0 1 9 6 5 1 9 7 0 1 9 7 5 1 9 8 0 1 9 8 5 1 9 9 0 1 9 9 5 2 0 0 0 2 0 0 5 2 0 1 0 2 0 1 5 0.6 0.8 T (a) 1 9 5 0 1 9 5 5 1 9 6 0 1 9 6 5 1 9 7 0 1 9 7 5 1 9 8 0 1 9 8 5 1 9 9 0 1 9 9 5 2 0 0 0 2 0 0 5 2 0 1 0 2 0 1 5 0.6 0.8 Q (b) 1 9 5 0 1 9 5 5 1 9 6 0 1 9 6 5 1 9 7 0 1 9 7 5 1 9 8 0 1 9 8 5 1 9 9 0 1 9 9 5 2 0 0 0 2 0 0 5 2 0 1 0 2 0 1 5 tween climate variability in distinct regions, which would be reflected by higher absolute correlation values. Specifically, in this work, we use a fixed link density of 0.5% in all window-specific climate networks and thus cannot make any statements about the overall strength of correlations. However, following previous results by Radebach et al. (2013), we may actually expect that the correlation threshold q |c|,0.995 used for establishing network connectivity in this work exhibits maxima whenever T shows a peak, thereby supporting the hypothesis of El Niño and La Niña episodes synchronizing global SAT variability by establishing teleconnections.
Regarding the latter observation, it is remarkable that previous works by other authors rather reported a reduction of connectivity associated with a breakdown of synchronization due to the large-scale climate disruption triggered by El Niño events . In fact, this observation has been recently used to develop network-based forecasting strategies for El Niño (Ludescher et al., 2013(Ludescher et al., , 2014. However, the apparent contradiction between the latter results and our observations can be resolved when taking the different approaches of network construction used in the respective works into account. Beyond their overall large-scale similarity, the temporal variability profiles of transitivity and modularity also exhibit some important differences. In particular, the strong 1982/83 El Niño is represented as a single long episode of reduced modularity values while being split into two rather distinct peaks in transitivity (see Fig. 2a and 2b). Given the known seasonal profile of El Niño peaking around Christmas, it is remarkable that the ONI index has remained high during a quite long period of time, indicating a single extended event even despite its temporary decay captured by T , but not quite by Q. This underlines that both measures actually capture different aspects of network organization that provide complementary information.
Another notable observation is related with the abrupt shift from El Niño to La Niña conditions in summer 1998, leading to a very fast reorganization of the global SAT field. The latter transition is reflected by some negative anomaly of T in summer 1998 in comparison with the "normal" background values of this measure, which presents a unique feature in the time evolution of network transitivity over the last decades that is not accompanied by any corresponding anomaly in Q.
Taken together, modularity and transitivity evolve similarly at larger time scales, but provide complementary viewpoints. High transitivity commonly coincides with the temporary appearance of densely connected structures in the functional climate network, which are typically well localized in space (Radebach et al., 2013). In turn, modularity captures the global connectivity pattern rather than primarily local features. Specifically, a low modularity value actually highlights more global connections in the climate network.
While network transitivity and modularity present two key topological network characteristics, functional climate networks are systems embedded in geographical space. Thus, the spatial placement of nodes and links (which is disregarded by topological characteristics) can play a pivotal role in network structure formation (Radebach et al., 2013). In order to address this aspect, we present the temporal evolution of the global average link distance d in Fig. 2c. Notably, this measure exhibits more irregular variability with a less clear distinction between "background level" and "anomalies" associated with different types of climate disruptions than the previously studied two topological characteristics. Yet, the general behaviour of d resembles that of network transitivity in the sense that ENSO-related peaks often co-occur in both measures. This indicates that strong El Niño and La Niña episodes do not exclusively trigger shortrange (localized) connectivity (high T ), but also global teleconnectivity (high d ), which is in line with contemporary knowledge on the large-scale impacts of both types of ENSO phases. Notably, this result is in agreement with previous qualitative results of Radebach et al. (2013) on the link distance distribution of global SAT networks.
From the results discussed above, we tentatively conclude that in order to distinguish globally influential ENSO events from episodes of minor (or more regional) relevance, a combination of modularity and average link distance can be useful, taking a holistic view in studying the differential imprints of different types of ENSO phases. We will recall this strategy when discussing the effects of volcanic eruptions on network organization at a global scale.

Spatial patterns of network connectivity
The above analysis of global network properties has largely confirmed some known effects of certain ENSO phases on the spatial co-variability structure of the global SAT field. Drawing upon the insight that topological and spatial network measures can provide different perspectives on the corresponding network patterns, we now turn to investigating the geographical characteristics of the generated functional climate networks. Specifically, following recent observations that climate network properties distinguish between the EP and CP flavours of both, El Niño and La Niña (Wiedermann et al., 2016), we are interested in the question how the associated (tele-) connectivity structures are manifested in the respective spatial fields of degree and average link distance. For this purpose, Fig. 3 shows composite plots of the spatial patterns exhibited by both network properties during the different types of ENSO phases, thereby averaging the local network properties over all time windows that are classified as showing either of these situations (see Tab. 2).
The left panels of Fig. 3 display the respective mean degree fields for the different types of ENSO periods. As ex-pected, we observe a particularly strong deviation from a homogeneous pattern during EP El Niños (Fig. 3a), while the degrees in the eastern-to-central tropical Pacific are only slightly larger than in the rest of the network during time windows without El Niño or La Niña conditions (Fig. 3i). This general behaviour is expected from previous studies (Wiedermann et al., 2016). Still, the observed degree patterns alone do not allow us to distinguish between a local or global phenomenon. For this purpose, the right panels of Fig. 3 show the corresponding mean average link distance fields for each type of situation. Elevated values of this measure in the typical ENSO region are present in case of all four possible types of episodes, indicating that both flavours of El Niño and La Niña actually generate additional connections in the tropical Pacific that span relatively large distances.
Analyzing the composite maps of the average link distance in more detail, it is important to note that beyond the ENSO region itself, additional parts of the globe exhibit elevated values. This indicates the presence of localized teleconnections that possibly link climate variability in the latter regions with ENSO. Specifically, EP El Niños (Fig. 3b) exhibit such teleconnections with Indonesia and Western Africa, which are also recovered for EP La Niñas (Fig. 3f). For CP El Niños (Fig. 3d), the d i field highlights a weak connection with Western Africa, but none with Indonesia. Similar but still weaker teleconnections can be observed for CP La Niñas (Fig. 3h).
Among the aforementioned patterns, the apparent teleconnection with Indonesia present during EP events but not during their CP counterparts is particularly interesting, as it is localized in the westernmost tropical Pacific. Thus, it connects eastern and western Pacific while not leading to marked long-distance connections in the central Pacific close to the dateline. One appealing explanation of this finding could be that the corresponding link is mediated via the Walker circulation (Ashok and Yamagata, 2009) and, thus, via airflow in higher altitudes rather than near-surface atmospheric circulation. However, it has to be noted that our analysis is based on cross-correlations only. The values can be severely affected by distinct temporal persistence properties of SAT in the eastern and western tropical Pacific, as pointed out by recent studies making use of modern causal inference methods (Balasis et al., 2013;Runge et al., 2014). Accounting for this effect in terms of replacing the correlation values by associated significance levels in the network generation step (Paluš et al., 2011) could provide a useful yet computationally demanding avenue for future research on this topic. From an impact perspective, the teleconnection suggested by our results is compatible with the documented increased likelihood of droughts in Indonesia during El Niño events (Diaz et al., 2001).
The apparent teleconnection with Western Africa spans a rather large spatial distance (about one third of the globe). In this context, Joly and Voldoire (2009) noted that "a significant part of the West African monsoon (WAM) interan- nual variability can be explained by the remote influence of El Niño-Southern Oscillation (ENSO)." This previously reported teleconnection could be responsible for the elevated average link distance over Western Africa especially during EP El Niños.
In general, climate variability in tropical regions is typically more likely to exhibit strong correlations than between tropics and extratropics, which is mainly due to the cellular structure of meridional atmospheric circulation that is effectively decoupling tropics and extratropics. In this regard, the omnipresent slightly elevated average link distance values in the polar regions are most likely data artifacts not corrected by our remapping procedure rather than indications of actual teleconnections.

Regionalized network characteristics
Global and local climate network properties as discussed above still provide only incomplete information on the effects of climate variability in different parts of the ENSO region on global SAT. To obtain further insights into this aspect, we now turn to analyzing the regionalized field measures introduced in Sect. 3.3 and study the specific connectivity associated with the Nino3.4, Nino3 and Nino4 regions in terms of degree and average link distance.
The corresponding results are summarized in Fig. 4. We observe that the relative magnitude of variations of regionalized degree and average link distance is even stronger than that of the global network properties transitivity, modularity and global average link distance discussed above. All measures exhibit episodes of very small values as opposed to such with much larger values, the latter often coinciding with El Niño and La Niña phases. Since the corresponding regions have been previously chosen for defining ENSO-specific indices, this result has been expected. Most importantly, degree and average link distance based characteristics exhibit strong positive correlations. Notably, for climatic events with predominantly local structure, we would expect a strong increase of k i but only weaker increase of d i in the region under study. Hence, our corresponding observations underline that ENSO-related climate impacts are not confined to the vicinity of the ENSO region, but are controlled by large-scale teleconnections.
Since the different ENSO regions show partial overlap (cf. Fig. 1), the results obtained for the individual regions exhibit a high degree of similarity. However, regarding specific El Niño or La Niña episodes, comparing the corresponding signatures for the Nino3 and Nino4 regions still allows attributing these events to more Eastern Pacific or Central Pacific types. For example, the strong 1997/98 El Niño is reflected by very high values of the regionalized degree for the Nino3 and Nino3.4 regions, but relatively weak signatures in the more western Nino4 region, which is consistent with its classification as an EP type event.
Examining the time evolution of all six regionalized network measures in some detail, it is notable that between 1978 and 1982, there has been considerable variability in all measures pointing towards an episodic presence of teleconnections even though none of the time windows was classified as an El Niño or La Niña episode according to the ONI. Moreover, we find that before the year 2000, clear peaks can always be observed in all properties as alternating with periods of low values. In turn, during the last about 15 years, we rather find strong variability without any low background level, with peaks occurring almost annually, with the exception of 2013 and 2014. This change in the overall temporal variability pattern of our regionalized network measures might point to some fundamental changes in the spatiotemporal organization of global SAT, either due to some not yet identified mode of natural variability or as a result of external interference. An attribution of this observation is, however, beyond the scope of the present work.

Volcanic eruptions
Besides ENSO variability, strong volcanic eruptions have been identified as causes of marked disruptions in climate network properties in earlier studies (Radebach et al., 2013). In this context, the application of the complementary viewpoints as used in this work for further characterizing the impacts of such eruptions promises interesting additional insights.
Regarding the global network properties, let us turn back to Fig. 2. As already emphasized in our discussion on the corresponding imprints of different ENSO phases, anomalies in transitivity and modularity need to be interpreted differently in terms of global versus more regional changes in climate network connectivity. While EP El Niños and La Niñas (but not their CP counterparts) consistently show peaks in transitivity coinciding with breakdowns in modularity, a similar signature has been found in the aftermath of the Mount Pinatubo eruption, suggesting that this event has affected the climate system globally. However, when comparing these topological network characteristics with the spatial network property of global average link distance d , we find a marked difference. Specifically, for ENSO-related climate disruptions, both d and T show the same direction of changes (i.e., peaks) with only few exceptions. In turn, we observe an opposite behaviour of both measures in 1993, which corresponds to the time windows where the cooling effects following the Mount Pinatubo eruption should have taken their maximum (McCormick et al., 1995). Hence, unlike for ENSO-related disruptions, the peak in transitivity and simultaneous drop in d indicate that the effects of the volcanic eruptions have rather been regionally confined. The latter is consistent with the hypothesis of elevated correlations in the region that has been most directly affected by the associated cooling trend following the eruption. Based on this observation, we suggest that using the global average link very strong strong moderate weak weak moderate strong very strong in conjunction with network transitivity and modularity enables us to discern disruptive events with global effects (strong ENSO phases) from those exhibiting more regional impacts (volcanic eruptions).
In general, one notable difference in comparison with ENSO-related impacts is that large-scale effects of volcanic eruptions on global SAT teleconnectivity can be observed only after a sufficiently large amount of aerosols have entered the stratosphere (Robock, 2000). Accompanying the resulting time shift between trigger event and response, we may also need to consider a spatial shift of the most affected region as compared to the location of the volcano. In the following, we apply our regionalization procedure described in Sect. 3.3.2 to studying the impacts of the Mount Pinatubo, Mount Agung and El Chichon eruptions. In order to avoid interference with the effects of ENSO events, the ENSO-big region depicted in Fig. 1 are excluded from the corresponding computations. The results obtained from this analysis are summarized in Fig. 5.
The largest of the three considered eruptions (Mount Pinatubo) had a global cooling effect and has left clearly visible signatures in all considered global network measures as discussed above. Some months after the eruption, a large region of elevated network connectivity has been established, which covers essentially all of the western tropical Pacific (Fig. 5c). The temporal evolution of the average degree in the region around Mount Pinatubo displays an abrupt rise about half a year after the eruption, then a constantly high value for about one year (the common residence time of volcanic aerosols in the stratosphere) before dropping again back to its previous level (Fig. 5a). The region with the highest degrees is shifted northward with respect to the location of the volcano (Fig. 5c). When computing the average degree for this region, we observe an even stronger rise of the regionalized degree than for the region surrounding the volcano (Fig. 5b).
The Mount Agung eruption exhibits similar, but weaker, patterns in the respective region (Fig. 5f). However, the region with the highest degree is shifted south-westward. The average degree in the region surrounding Mount Agung only shows weak changes after the eruption (Fig. 5f), as opposed to a somewhat sharper increase in the shifted region, with the peak effect occurring significantly faster after the beginning of the eruption than in case of the Mount Pinatubo eruption (Fig. 5e). However, it should be noted that we can already observe the beginning of some upward trend in regionalized degree before the actual eruption, pointing to a possible interference with normal natural variability.
Unlike the two other volcanic eruptions, the degree field in the period succeeding the El Chichon eruption showed hardly any marked changes (Fig. 5i). Consequently, we also do not observe any marked signature in the temporal variability profile of the regionalized degree in the surrounding of the volcano (Fig. 5g). Instead of a peak shortly after the eruption, we actually find a clear drop of the average degree. However, given that El Chichon is located relatively close to the ex-tended ENSO region, we cannot rule out that this could be an effect of the strong El Niño event occurring shortly after the eruption and eventually even being partially triggered by the latter (Khodri et al., 2017). In general, previous studies indicate that the El Chichon eruption caused a much weaker summer cooling than the Mount Agung eruption (Man et al., 2014), which could also explain its absent signature in our analysis.

Conclusions
We have used functional climate networks constructed from spatial correlations of daily global surface air temperature (SAT) anomalies to analyze the global impact and teleconnections of past El Niño and La Niña events as well as volcanic eruptions. By making use of the global network property of modularity, we have found that at least the East Pacific flavours of such events lead to a global reconfiguration of SAT variations. Considering the global average link distance as a complementary spatial network characteristic, we have identified distinct qualitative differences between the imprints of these ENSO periods and the Mount Pinatubo eruption in global SAT patterns.
Using composites of the degree and average link distance fields, we have identified hallmarks of distinct ENSO teleconnections in the climate network structure, especially such linking the eastern tropical Pacific with Indonesia and West Africa during East Pacific El Niños, both of which have already been reported in previous studies (Diaz et al., 2001;Joly and Voldoire, 2009). By making use of a regionalization procedure applied to these two fields of local network properties, we have introduced a simple yet effective tool to unveil the differential roles of different regions in the tropical Pacific in establishing teleconnections during different El Niño and La Niña events.
Finally, we have analyzed the global and local connectivity properties of SAT-based climate networks in the aftermath of the strongest recent volcanic eruptions of Mount Pinatubo, Mount Agung and El Chichon. In particular, while the Mount Pinatubo eruption has been confirmed to exhibit marked impacts on global SAT, its dominating effect was rather regional (i.e., it did not trigger long-range teleconnections detectable by our approach).
While most of the results presented in this work rely on a qualitative analysis of temporal changes in the climate network properties, additional statistical quantification of their relationship with existing indicators of ENSO variability and volcanic eruptions' impacts might further strengthen our findings. Regarding ENSO, many previous studies have attempted to utilize the spatial patterns of SST anomalies to define corresponding index variables. However, the corresponding classifications of El Niño and La Niña phases reached only partial consensus, which was in fact the motivation of the work of Wiedermann et al. (2016) presenting Figure 5. Time series of regional mean degree and associated degree field (excluding the ENSO-big region indicated by white color) for the three main volcanic eruptions during the study period: (a-c) Mount Pinatubo, (d-f) Mount Agung, (g-i) El Chichon. In the degree maps shown in panels (c), (f) and (i), blue dots mark grid points within a radius of 5 • around each volcano, which have been used to define the regionalized degrees shown in panels (a), (d) and (g), respectively. Red dots indicate spatially shifted regions of the same size where the largest changes to the degree field have been observed. These regions serve as the basis for computing the regionalized degrees shown in panels (b), (e) and (h), respectively. Purple vertical lines indicate the timing of the respective eruptions, whereas green vertical lines indicate the midpoints of the time windows exhibiting the strongest signature in the regionalized network properties. The time series have been restricted to ±10 years around the date of the respective eruption. Background colours indicate the corresponding ENSO strength as in Fig. 2. climate network transitivity as a useful and consistent index. Going one step further, one might easily quantify, for example, the correlation between transitivity and other (global or regionalized local) network characteristics. However, in our opinion, the particular value of the present work lies in identifying specific properties that are not exhibited by the former (as well as other not network-related indices). In this context, there is no established benchmark that could be used for further testing the significance of our results.
In turn, regarding the effects of volcanic eruptions, the respective regionalized degrees for the spatially shifted "major impact regions" of both, Mount Pinatubo and Mount Agung, exhibited their overall maximum values among all time windows studied in this work in the aftermath of the associated eruptions. This indicates a very high significance of our corresponding results. Note that, however, we did not succeed in finding any comparatively strong impact signature in the climate network properties after the eruption of El Chichon, as well as after other strong volcanic eruptions of the past about 70 years (not shown). We relate the latter finding to the generally lower magnitude of the respective perturbations (in terms of a smaller amount of climate-active volcanic aerosols injected into the stratosphere). Moreover, some of the other major eruptions (e.g., the Mount St. Helens eruption in 1980) appeared in the extratropics rather than the tropics. Together with the different seasonality of these events, this could imply different effects on regional and global temperature patterns, similar as shown recently for global monsoon precipitation (Liu et al., 2016).
In summary, our study confirms that ENSO does not only have a strong local effect on SAT in terms of coherent SAT trends in the tropical Pacific associated with a spatially confined increase of network connectivity (Radebach et al., 2013), but also dynamically reconfigures climate variability globally by triggering teleconnections especially with other tropical regions. In this regard, one possible mechanism could involve the modulation of monsoons by strong El Niño and/or La Niña periods, which could be further modulated by volcanic eruptions (Maraun and Kurths, 2005). Confirming this hypothesis in the context of climate network studies would, however, require much more elaborated approaches than those used in the present work, and is therefore outlined as a subject of future research.

Appendix A: Comparison of different modularity estimation algorithms
In Fig. A1, we show the results of five algorithms to estimate the community structure of our functional climate networks in terms of the resulting modularity values: fast greedy (Clauset et al., 2004), infomap (Rosvall and Bergstrom, 2008), label propagation (Raghavan et al., 2007), leading eigenvector (Newman, 2006a) and WalkTrap (Pons and Latapy, 2006). Further details motivating the choice of Walk-Trap as a reference algorithm in the body of this paper are provided in the figure caption.
Data availability. All data used in this work are publicly available from NCEP/NCAR. In addition to the SAT data, we have made use of the monthly Ocean Niño Index (ONI) values as provided by NCEP (2017).
Code availability. All code used in this work has been written in Python and published under GPLv3 license as a GitHub repository (Kittel, 2017). Detailed information for the reproduction of the results of this paper can be found there. While the published code was originally designed to produce these specific result, we kept it rather general with further extensions in mind. Thereby, it can be used as a starting point for future evolving network research, as it provides some basic structures that are needed for evolving network analysis, for example an interface for hdf5 (a high-performant file type for data storage) and automatic parallelization using mpi.
Author contributions. TK, CC, NL, FR and RVD designed the analysis. TK, CC and NL conducted the analysis. TK, CC and RVD prepared the manuscript. TP, FR, JK and RVD supervised the analysis and revised the manuscript and the interpretation of the obtained results.
Competing interests. The authors declare no conflict of interest.
1 9 5 0 1 9 5 5 1 9 6 0 1 9 6 5 1 9 7 0 1 9 7 5 1 9 8 0 1 9 8 5 1 9 9 0 1 9 9 5 2 0 0 0  Figure A1. Comparison of estimated modularity values for the functional climate networks obtained for running windows as described in the main text. We use five different algorithms for detecting the underlying community structure. Since modularity estimation resorts to a numerical maximization problem, higher values indicate better results. Visual comparison reveals that the leading eigenvector and Walk-Trap algorithms outperform the others regarding this criterion. Since the leading eigenvector algorithm suffers from intermittent modularity breakdowns, possibly indicating numerical instabilities, we use the WalkTrap method in this paper.