Evolving climate network perspectives on global surface air temperature effects of ENSO and strong volcanic eruptions

Abstract

Episodically occurring internal (climatic) and external (non-climatic) disruptions of normal climate variability are known to both affect spatio-temporal patterns of global surface air temperatures (SAT) at time-scales between multiple weeks and several years. The magnitude and spatial manifestation of the corresponding effects depend strongly on the specific type of perturbation and may range from weak spatially coherent yet regionally confined trends to a global reorganization of co-variability due to the excitation or inhibition of certain large-scale teleconnectivity patterns. Here, we employ functional climate network analysis to distinguish qualitatively the global climate responses to different phases of the El Niño–Southern Oscillation (ENSO) from those to the three largest volcanic eruptions since the mid-20th century as the two most prominent types of recurrent climate disruptions. Our results confirm that strong ENSO episodes can cause a temporary breakdown of the normal hierarchical organization of the global SAT field, which is characterized by the simultaneous emergence of consistent regional temperature trends and strong teleconnections. By contrast, the most recent strong volcanic eruptions exhibited primarily regional effects rather than triggering additional long-range teleconnections that would not have been present otherwise. By relying on several complementary network characteristics, our results contribute to a better understanding of climate network properties by differentiating between climate variability reorganization mechanisms associated with internal variability versus such triggered by non-climatic abrupt and localized perturbations.

1 Introduction

The empirical analysis of climate data is fundamental to understand the evolution (and subsequently 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, which have been observed at distinct locations distributed around the globe. Their common properties include long-range spatial and often also temporal correlations [1], interactions at and among multiple scales [2] and nonlinearity [3]. 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 (EOF) or maximum covariance analysis [4], it has recently been realized that these methods exhibit fundamental intrinsic limitations, including their linearity and associated condition of pairwise orthogonal patterns (e.g. [5]). As a consequence, the traditional view that the corresponding decompositions of global spatio-temporal co-variability patterns actually provide dynamical (or at least statistical) modes that unambiguously coincide with specific key climatic processes has been challenged [6]. Along with modern methodological developments in complex system theory and computer sciences, there is a growing body of literature demonstrating examples of successful applications of nonlinear methods to climate problems and their ability to unveil relevant dynamical characteristics that may be hidden to traditional analysis techniques [7,8,9,10].

One specific example of such methodological developments are complex network representations of climate variability [11,12,13,14,15,16,17,18,19]. Although they share some common roots with classical techniques like EOF analysis, they generalize the corresponding scope and can potentially relieve some of the aforementioned concerns [20]. In this work, we focus on the so-called functional climate network analysis, in which the individual grid points or cells are considered as nodes of a spatially embedded graph. Connections among these nodes are established according to dynamical similarities between the individual (local) climate time series [11, 13, 21]. By construction, the network structures thus obtained highlight essential statistical interrelationships among spatio-temporal climate data [13].

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 [12]. Moreover, climate networks have been used to identify dipole patterns which represent pressure anomalies of opposite polarity appearing in two different regions simultaneously [22]. The study of the coupling structure between interdependent climate variables [23], the temporal evolution and teleconnections of the North Atlantic Oscillation [24, 25], the distinction of different types of El Niño phases [26, 27], and even the prediction of El Niño occurrences [28, 29] and magnitudes [30] have also been subjects of corresponding recent investigations. Many of the aforementioned methodological achievements have been integrated in open source software packages [31, 32] contributing to the increasing use of functional network analysis in climatological studies [21].

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 [14, 33, 34]. Here, a community is a subset of densely connected nodes which exhibit only few interactions with the rest of the network [35, 36]. In a climate network context, communities would ideally have some climatological interpretation. Specifically, Tsonis et al. [14] 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 [36] can also be used to detect multi-variable clusters [37]. Conceptually related to the aforementioned line of research are approaches attempting to define climate areas, i.e., regions with coherent climate variability, based on spatial clustering of grid cells in a fully connected weighted climate network [38] or other multivariate analysis techniques [10, 39].

In this paper, we present an attempt to further clarify the distinct nature of modifications in the climate network configuration arising along with either different types of El Niño and La Niño episodes or strong volcanic eruptions. Previous work has demonstrated that global surface air temperature (SAT) anomaly based networks exhibit clear temporary changes in some global network characteristics such as transitivity, average path length or global clustering coefficient [26] for the East Pacific types of El Niño and La Niño [27] while not for their Central Pacific counterparts, but have also shown similar anomalies in the aftermath of the Mount Pinatubo eruption. To further address the aforementioned aspects in greater detail and differentiate more clearly between different mechanisms affecting the connectivity reorganization processes in evolving climate networks, we complement the previously employed network transitivity metric with information on the modular organization of the climate network and the presence of long-range teleconnections. Specifically, we study the co-variability structure of the global SAT field for running windows in time in two different ways: On one hand, we consider spatial fields of two simple network properties that represent the number of strong statistical connections associated with a given grid point, as well as the average spatial distance between the connected grid points (i.e., an indicator for possible teleconnectivity). On the other hand, the associated temporal variations in global network connectivity properties are traced by a suite of scalar-valued characteristics quantifying (i) the transitivity of strong correlations, (ii) the modular organization of the network, and (iii) the average distance between mutually connected grid points. We show that these three measures capture complementary aspects of the spatial organization of surface air temperature co-variability and when combined allow distinguishing the effects of large-scale climate disruptions into regional and teleconnectivity effects. Specifically, we provide empirical evidence that strong ENSO episodes are able to generate both enhanced regional as well as long-distance connectivity patterns simultaneously, whereas the strongest volcanic eruptions since the mid-20th century rather had predominantly regional effects. Moreover, the conditions under which specific ENSO phases result in qualitative large-scale connectivity reorganization of the SAT network are further detailed by analyzing the network patterns in dependence on both ENSO diversity (i.e. East Pacific versus Central Pacific flavors) and transitional complexity (i.e. the succession of ENSO phases).

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. Additional information on the robustness of our findings along with modified analysis settings is provided in two brief appendices.

2 Climatological background

2.1 El Niño–Southern Oscillation

Fig. 1

Main regions of interest used within this paper. Sets of blue dots labeled with “El Chichon”, “Agung” and “Pinatubo” indicate grid points within a 5$^\circ$ radius around the corresponding volcanoes. The numbers 3, 3.4 and 4 identify the corresponding Nino regions (cf. Table 1) commonly used for defining characteristic indices of ENSO variability (marked with green background color, red line, and blue background, respectively). The region “ENSO-big” (dashed green line) will be removed from the complete global data set when analyzing the spatial imprints of volcanic eruptions to ensure that ENSO-related regional effects are excluded

Among the dominant spatio-temporal variability 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 [40, 41]. During the positive phase (El Niño) of this complex oscillation of the coupled atmosphere–ocean 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 anomaly results in marked shifts of key atmospheric pressure systems [42]. Such shifted pressure patterns cause large-scale changes in the atmospheric circulation, leading to prominent shifts of precipitation patterns. More specifically, tropical Pacific SST anomalies associated with different ENSO phases alter the location of main tropical convection zones, which modify the location and intensity of diabatic heating as the main driver of tropical and extratropical circulation patterns. In this context, the response of atmospheric circulation to El Niño or La Niña events is not restricted to the tropical Pacific, but also affects the extratropics via the excitation of atmospheric waves migrating towards the extratropics and/or interacting with midlatitude dynamics [43, 44]. It has been shown that the resulting effects of both ENSO phases can be observed in remote regions including North and South America, Africa, the Indian subcontinent, and even Antarctica [45,46,47,48,49,50].

The interannual variability of ENSO is dominated by some irregular oscillations with a period of 2–7 years and remarkable variations in the associated characteristic frequencies and amplitudes of both temperature and pressure anomalies [40]. 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 [51] (see Fig. 1 and Table 1). In this work, we will utilize the so-called Oceanic Niño Index (ONI) for differentiating between different phases of ENSO. It is defined as the running three-month mean SST anomaly for the Niño 3.4 region ($5^{\circ }\hbox {N}$–$5^{\circ }\hbox {S}$, $120^{\circ }$–$170^{\circ }\hbox {W}$) in comparison with centered 30-year base periods that are updated every 5 years [52]. When the ONI exceeds $0.5 ^{\circ }\hbox {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 ^{\circ }\hbox {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 flavors (see [27] and references therein). The first type is the canonical or East Pacific (EP) El Niño [53, 54], the main SST anomaly pattern of which is located 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 [42]. Notably, both spatial structures (EP and CP) can also be observed in the context of La Niña, yet commonly appearing less distinct in the corresponding SST anomaly fields. Noticing that there had been contradictory classifications in the literature for some past ENSO events, Wiedermann et al. [27] recently presented a new indicator for the ENSO flavor based on functional climate networks (i.e., the same methodology that is also used in the present work). In the remainder of this paper, we will follow their classification for consistency reasons, which is summarized in Table 2.

2.2 Volcanic eruptions

Besides distinct ENSO episodes and their known global climate impacts, another type of events that can substantially affect climate at considerable 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 temperature trends that may last up to several years depending on the magnitude of the event [55]. During strong volcanic eruptions, large amounts of sulphate aerosols are injected into the stratosphere, which have a typical residence time of the order of one year. The presence of the volcanic aerosol cloud has differential effects on temperatures in the covered region, including tropospheric cooling (via a reduction of solar radiation reaching the surface) and stratospheric warming (via an enhanced absorption of radiation). These modifications in the radiation balance and atmospheric chemistry may affect climate variability in different ways, depending on the spatial location (specifically, its geographical latitude) of the volcano and the season of the eruption. Details on the corresponding processes can be found in [55]. As a consequence, volcanic eruptions can again cause changes of precipitation and temperature patterns from synoptic (weather) time scales to relatively persistent multi-annual effects [55] and even trigger long-lasting climate disruptions [56, 57] and corresponding societal impacts [58].

Table 1 Overview on different regions commonly used for defining characteristic temperature-based indices associated with ENSO variability

It is interesting to note that the respective climate disruptions caused by ENSO and volcanic eruptions may not be mutually independent. In fact, there exists some evidence for volcanic eruptions having different effects on El Niño depending on the latitude and timing of the eruption [59, 60] or even serving as triggers of El Niño events [61], which yet seems to be poorly represented in existing tropical Pacific paleoclimate proxies [62]. The specific climatological questions related with volcanic effects on ENSO and, more generally, climate variability, including the episodic synchronization between ENSO and Indian monsoon rainfall in the presence of strong volcanic activity [63], have recently triggered enormous research efforts among the climate modeling community.

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 entering 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 [64], took place between April and September 1991 in the Philippines, followed by the Mount Agung eruption in Indonesia (February 1963 to January 1964) and the El Chichon eruption in Mexico (March to September 1982) (see Fig. 1).

Table 2 Classification of past ENSO episodes into the four types Eastern Pacific (EP) El Niño, Central Pacific (CP) El Niño, EP La Niña, and CP La Niña as proposed by [27].

3 Data and methods

3.1 Data

We use daily mean surface air temperature (SAT) data (at sigma level $\sigma =0.995$) from the National Center for Environmental Prediction (NCEP) and National Center for Atmospheric Research (NCAR) Reanalysis I project [65, 66]: The data cover the years 1948–2015 at a global grid with equi-angular spatial resolution of $2.5^{\circ }$ in both latitude and longitude, thus comprising 10, 512 individual temperature time series. 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 close to the poles than around the equator. This 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 [11, 67, 68]. As an alternative, we follow here the approach of Radebach et al. [26], where the original data have been remapped onto a grid with a much higher spatial homogeneity. Specifically, we use an icosahedral grid as described in [69], 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.

Fig. 2

Time series of climate network a transitivity ${\mathcal {T}}$, b modularity ${\mathcal {Q}}$ and c global average link distance $\left\langle \left\langle d\right\rangle \right\rangle$. Background colors highlight different ENSO phases (red: El Niño (EN), blue: La Niña (LN)) according to the Oceanic Niño Index (ONI), with opacity representing the corresponding index value. Ticks on the time axis indicate the 1st of January of a given year; all values are shown according to the midpoint dates of the respective time windows

In [26], the time series associated with each new grid point have been determined based upon the values at the respective four surrounding grid points forming a spherical rectangular cell of 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 neighbors 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 precipitation) and its resulting spatial smoothness, we do not expect the time series resulting from both algorithmic variants to differ markedly, which is also supported by the results on the time-dependent network transitivity measure (see below) in [26] in comparison with those obtained in our study (see Sect. 4.1.1, Fig. 2a).

Notably, differences between the two interpolation strategies can be expected to arise mainly in regions with very asymmetric original grid cells, i.e., close to the polar where the spherical rectangular cells of the original data grid are extremely stretched. When using the four nearest grid points for interpolation in the current work, we may accept that this can lead to using four points from the same latitudinal circle, thereby overemphasizing zonal correlations while possibly discarding meridional ones. Despite the existence of physical processes linking tropical and polar climate variability (typically via stratospheric pathways implying a considerable time lag), we may argue that the strength of instantaneous co-variability in the tropics (being the main source of climate disruptions studied in this work) and the Arctic and Antarctic is very likely relatively small, so that the polar regions potentially affected by the aforementioned zonal correlation bias can be expected to not contribute much to our analysis results. A more detailed quantitative evaluation of this assumption is, however, beyond the scope of the present work. In terms of correlation based functional climate network properties studied in this work, a possible zonal correlation bias could be identified by employing directional linkage statistics, which for the case of complex networks embedded in a spherical geometry have recently been introduced and discussed by Wolf et al. [70].

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. To account for this problem of temporal co-occurrence between the effects of volcanic eruptions and ENSO events, there exist different possible strategies.

On one hand, we may attempt to remove the effect of ENSO variability from all considered SAT time series, commonly by considering the residuals of a corresponding linear regression [71, 72]. In this case, all further analyses would be based on pairwise partial correlations between SAT time series conditional on ENSO. Similar approaches of partial correlation or partial regression maps have been considered in previous climatological studies. To this end, we, however, have to recall that this approach would involve strong assumptions and some associated conceptual disadvantages. First, it is not clear which specific index of ENSO variability should be used for regression, and the results may depend on a corresponding choice. Second, using linear regression makes assumptions on the way ENSO may affect local SAT variability, which appear not quite justified for a complex nonlinear phenomenon like ENSO. Using other (nonlinear) functional parameterizations of those effects would not quite solve this conceptual problem while introducing additional numerical efforts. Third, ENSO itself is strongly entangled with variability at other time scales like annual or quasi-biennial ones [5, 73], so that a sufficient removal of its effects by simple regression appears unlikely. Finally, it may well be that the ENSO background state can pre-condition the response of the global climate system to specific volcanic eruptions [61, 63], which would not be accounted for in this approach.

On the other hand, given the aforementioned considerations on the partial regression method, we will approach the problem of separating the respective effects of ENSO and volcanic eruptions from a spatial perspective, considering that the ENSO phenomenon itself is (despite its global teleconnectivity) confined to a region which does not overlap spatially with any of the three major volcanic eruptions in the last decades. Thereby, removing all grid points from the greater ENSO region in the tropical Pacific Ocean allows us studying the spatiotemporal organization of SAT co-variability that is not explained by direct linkages with this region of the globe. Following this idea, 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 Table 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 [26].

3.2 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 [11, 15, 19, 21]. 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 associations 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.

Since there exists a growing body of literature on applications of functional climate network analysis, we refer the reader to some of the existing reviews [19, 21] for a more exhaustive description of the methodology. To this end, let us consider evolving climate networks constructed from sliding windows of daily SAT anomalies (with the local mean annual cycle being removed) with a width of $\varDelta d=365$ days and mutual offset $\delta$ between subsequent windows. In this work, we have generally set $\delta =1$ day.

For the resulting time series segments corresponding to a specific one-year time window, we compute the classical lag-zero pair-wise Pearson correlation coefficients $c_{ij}$ and identify the empirical $99.5\%$ quantile $q_{|c|,0.995}$ of all absolute values $|c_{ij}|$ for each time window (the effects of different thresholds will be briefly summarized in Appendix A). With this quantile, which is empirically determined individually for each time window and thus changes with time, we define the adjacency (connectivity) matrix ${\mathbf {A}}$ of the climate network for a window centered at time $d_{mid}$ as

$$\begin{aligned} A_{ij}(d_{mid}) = \varTheta (|c_{ij}|(d_{mid})-q_{|c|,0.995}(d_{mid}))-\delta _{ij}\,,\qquad \end{aligned}$$

(1)

where $\varTheta (\bullet )$ is the Heaviside function and $\delta _{ij}$ denotes Kronecker’s Delta. Notably, ${\mathbf {A}}\left( d_{mid}\right)$ contains the complete structural information on the evolving network representation of the global surface air temperature field.

Based on the thus obtained adjacency matrices, we consider the following network properties in dependence on time:

• The degree $k_i = \sum _{j=1}^N A_{ij}$ of a node i measures how densely this node is connected within the network. In case of a functional climate network, it provides a proxy for the importance (or centrality) of a certain grid point in the spatio-temporal correlation structure of the variable of interest. Network measures like the degree, which are defined individually for each node i, will be referred to as local network characteristics or spatial fields in the following.

• The local average link distance [21]

  $$\begin{aligned} d_i = \left\langle d_{ij}\right\rangle _{\left\{ j|A_{ij}=1\right\} }&= \frac{1}{k_i} \sum _{\left\{ j|A_{ij}=1\right\} } d_{ij}\nonumber \\&= \frac{1}{k_i} \sum _{j=1}^N A_{ij}d_{ij}\,, \end{aligned}$$

  (2)

  measures the mean spatial distance covered by all links associated with a given node i, where $d_{ij} = 2 D_{ij}/u_{Earth}$ with $D_{ij}$ being the geodesic distance between two nodes i and j and $u_{Earth}$ the circumference of the Earth. A low average link distance indicates very localized connections, while a high value points to a node with long-ranging spatial connections (teleconnections). Taking the average of $d_i$ over all nodes i of the network gives the global average link distance $\left\langle \left\langle d\right\rangle \right\rangle = \left\langle d_i\right\rangle _i$. For the sake of brevity, we will use the abbreviated term average link distance in this manuscript for describing both the local ($d_i$) and global ($\left\langle \left\langle d\right\rangle \right\rangle$ versions of this property as long as the corresponding context is unambiguous.

• The transitivity

  $$\begin{aligned} {\mathcal {T}} = \frac{\sum _{i,j,k=1}^N A_{ij} A_{jk} A_{ki}}{\sum _{i,j,k=1}^N A_{ij} A_{jk}} \end{aligned}$$

  (3)

  quantifies how strongly the connectivity of a network is clustered, i.e., the fraction of cases in which the presence of two links between nodes i and j as well as i and k also implies a link between j and k [26, 74]. Like the global average link distance $\left\langle \left\langle d\right\rangle \right\rangle$ (but unlike degree and local average link distance), ${\mathcal {T}}$ does not define a field, but returns one single scalar value for each network.

• The modularity [75]

  $$\begin{aligned} {\mathcal {Q}} = \frac{1}{2m} \sum _{ij}\left( A_{ij} -\frac{k_i k_j}{2m}\right) \varDelta _{ij} \end{aligned}$$

  (4)

  (where m is the total number of links and $\varDelta _{ij}$ an indicator function informing whether or not two nodes i and j belong to the same subgraph in a given partition of the network) measures 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 context of a climate network, modularity provides a single scalar-valued 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 individual subgraphs in the partition that maximizes the modularity ${\mathcal {Q}}$ are called communities. The higher the modularity, the more split-up (or modular) the network. There are a plethora of community detection algorithms [36], which provide different solutions and, hence, different resulting modularity values. Here, we employ the WalkTrap method [76], which has been found to exhibit comparably high values of modularity and relatively stable values in case of strongly overlapping windows (see Appendix B for details). We note that due to the relatively large computational efforts associated with the community detection, we have performed this analysis step not for all possible time windows of observations with mutual offset of $\delta =1$ day, but only for every 15-th window (i.e. with an effective offset of $\delta =15$ days between subsequent values).

Table 3 Overview on the three largest volcanic eruptions during our study period, along with the spatial location of the volcano and the center of functional climate network connectivity changes (degree anomalies) in the aftermath of the eruption (see Sect. 3.3.2 for details on how the latter have been obtained)

It is important to highlight that different from other studies, we do not impose the same threshold value to the absolute correlation coefficients for constructing our functional climate networks (which would allow the networks’ link densities to vary in time), but initially select the link density and let the threshold vary accordingly. This strategy of a fixed link density has the important advantage that the network characteristics obtained for different time windows are comparable, while they could otherwise exhibit an intrinsic dependence on the link density. This applies, for example, to the modularity – if we prescribe a certain network organization principle (e.g. a random graph or preferential attachment procedure), varying the link density can lead to quite different modularity values. Similarly, since most links in a climate network span short distances due to spatial autocorrelation effects of the global SAT field, elevating the link density (and thereby lowering the absolute correlation threshold) will give rise to more and more longer links to emerge, thereby systematically elevating the average link distance. To exclude such implicit density effects and focus on the rearrangement of connectivity between different spatial scales, all results reported in this work will be based on network configurations with the same fixed link density of 0.5% (unless stated otherwise). This setting is in agreement with previous studies [26, 27].

3.3 Regionalization of field measures

As already noted 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 ${\mathcal {X}} \subseteq \left\{ 1, \dots , N \right\}$ representing a certain part of the globe, a regionalized version of the degree would be given as

$$\begin{aligned} k_{{\mathcal {X}}} = \frac{1}{|{\mathcal {X}}|}\sum _{i \in {\mathcal {X}}} k_i\,, \end{aligned}$$

(5)

where $|{\mathcal {X}}|$ denotes the number of nodes in the considered set. Note that we split here the already constructed network representation into subsets of nodes and evaluate their respective features instead of constructing separate networks for each subset. 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 [23], since unlike $k_{{\mathcal {X}}}$, the latter exclude contributions due to links between nodes within ${\mathcal {X}}$ in their definition.

For the average link distance, the corresponding regionalized property $d_{{\mathcal {X}}}$ can be defined in full analogy.

Below, we detail some reasonable choices for ${\mathcal {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.

3.3.1 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 regions (Nino1+2, 3, 4 and 3.4) have been previously 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 applying Eq. (5). Thereby, we obtain a set of eight new scalar-valued characteristics: $k_{\text {Nino1+2}}$, $d_{\text {Nino1+2}}$, $k_{\text {Nino3}}$, $d_{\text {Nino3}}$, $k_{\text {Nino4}}$, $d_{\text {Nino4}}$, $k_{\text {Nino3.4}}$ and $d_{\text {Nino3.4}}$. 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.

3.3.2 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 structures in the global SAT field possibly induced by these eruptions, we need to aggregate 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^{\circ }$ around the location of each volcano (see Table 3) as basis for the regionalization procedure of $k_i$. This leads to the three observables $k_{\text {Pinatubo}}$, $k_{\text {Agung}}$ and $k_{\text {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 geographical neighborhood of a volcano 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 surface cooling to be expected. Accordingly, the induced anomalous (tele-) connectivity 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^\prime _{\text {Pinatubo}}$, $k^\prime _{\text {Agung}}$ and $k^\prime _{\text {Chichon}}$.

For determining the respective optimal spatial shifts for the main climate network responses to each of the three volcanic eruptions, we have examined the resulting degree fields for time windows with starting dates succeeding the individual eruptions by up to one year, seeking for the timing and position of the strongest anomalies in the degree field that could be attributed to each eruption. To understand this temporal offset, it should be noted that although a volcanic eruption may start relatively abruptly and generate the strongest injection of sulphate aerosols into the stratosphere during its initial stages, its larger-scale atmospheric effects due to the modified atmospheric chemistry and thereby affected radiation balance commonly become effective only with a considerable delay of several months up to two years after the start of the eruption [55, 64].

In addition to examining the resulting degree fields, the typical lower stratospheric wind fields during the period of each eruption have been considered as a consistency check for compatibility with a physically plausible response. Regarding possible latitudinal offsets, during the initiation of the different eruptions we may have expected no strong meridional wind component at about 15–20 km altitude for the Mount Agung and El Chichon eruptions, but some northward motion during the Mount Pinatubo eruption due to the prevailing seasonality of the stratospheric Brewer-Dobson circulation [77]. The zonal displacement of volcanic aerosols during tropical eruptions is being controlled by the quasi-biennial oscillation (QBO), with lower troposphere easterlies in the aftermath of the Mount Agung and El Chichon eruptions as opposed to weak westerlies during the Mount Pinatubo eruption [78, 79].

The results obtained from the described qualitative analyses provide consistent information on the expected as well as the empirically observed regions exhibiting the strongest effects of the three considered volcanic eruptions on the regional organization of network connectivity. In summary, our analysis is thus able to provide the approximate locations of the strongest degree anomalies in the SAT network (which do not necessarily coincide with the strongest surface cooling effects) as reported in Table 3, which are subsequently used for further study on the network connectivity reconfiguration during the aftermath of the three eruptions.

4 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 modified connectivity patterns induced by strong volcanic eruptions.

4.1 El Niño–Southern Oscillation

We start with investigating the 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. [26, 27, 80,81,82,83] 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.

4.1.1 Global topological network properties: transitivity and modularity

For almost the same analysis setup as used in our present study (same data set except for the last few years on record, same climate network link density of 0.5%, same icosahedral grid with only minor modification of the interpolation scheme), Radebach et al.  [26] demonstrated a striking co-variability between the threshold value of the absolute correlation coefficients used for the functional climate network construction, the network transitivity and the conceptually related global clustering coefficient, and the average path length, all of which exhibited exceptional values along with certain ENSO phases and the aftermath of the Mount Pinatubo eruption. We refer to Fig. 5 of the aforementioned reference for details and refrain from reproducing the corresponding results at this point. Most notably, correlation threshold and network transitivity were found to exhibit a very strong linear correlation of about 0.85 at zero and one month time lags, cf.  Fig. 7 in [26]. Following upon the latter result, we expect that the time dependence of the correlation threshold does not add much complementary information that is not provided by the corresponding behavior of the network transitivity.

Importantly, the network transitivity ${\mathcal {T}}$ has been shown by [27] to systematically discriminate between the EP and CP flavors of both El Niño and La Niña. While the mentioned reference used an area-weighted version of ${\mathcal {T}}$ on the original regular latitude–longitude grid of the reanalysis data and included information on the total pairwise correlation strength instead of just binary adjacency information, we follow here the earlier approach of [26] in using a remapping of the data onto an icosahedral grid along with unweighted network properties. 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. Note that Fig. 2 shows the results for different network measures in dependence on the window midpoint, while [27] used the endpoint, leading to a 6-month shift between the respective plots.

Radebach et al. [26] found that the emergence of transitivity peaks along with EP El Niño and La Niña phases is accompanied by a distinct reorganization of network connectivity, which involves the simultaneous appearance of (i) more short-range connections (leading to denser connectivity especially in the tropical Pacific, which reflects synchronous large-scale warming/cooling trends elevating the spatial autocorrelations in the region) along with (ii) more very long-distance links (describing emerging teleconnections). Necessarily, these two effects are compensated by a reduced number of links with intermediate distances because of the overall conservation of the number of links (fixed link density). Notably, the high transitivity appears to be mainly affected by (i), i.e. enhanced regional connectivity, because those short links are far more abundant than the long distance connections and therefore contribute more to the transitivity measure. In this regard, it is questionable if the network transitivity alone can differentiate between the two described effects.

To further quantify the strength of teleconnectivity in the global SAT field relative to local connections originating from spatial autocorrelations, we suggest that the network modularity ${\mathcal {Q}}$ can provide 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 large-scale organization of the global SAT network, i.e., a tendency towards more balanced co-variability in global temperatures across various spatial scales, particularly involving long-distance teleconnections.

Figure 2b indicates a pronounced negative correlation between transitivity and modularity for the investigated SAT networks. Specifically, almost all (see below for a more detailed discussion) time intervals that are characterized by elevated values of network transitivity actually exhibit a marked reduction in modularity, and vice versa. Consistent with previous findings [26], 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. Specifically, the drop in modularity appears compatible with the expected signature of emerging teleconnectivity that has been reported previously [26, 27] to occur along with EP events, but not CP phases—at least not to an extent detectable by network transitivity. Note that taken alone, this process would not necessarily imply a stronger synchronization between climate variability in distinct regions, which would be reflected by higher absolute correlation values. Specifically, in this work, we have consistently used a fixed link density of 0.5% in all window-specific climate networks and thus did not investigate the overall strength of correlations. However, following previous results [26] demonstrating a very strong positive correlation between absolute correlation threshold $q_{\left| c\right| ,0.995}$ used for establishing network connectivity and network transitivity at the studied link density, we may actually expect that the threshold exhibits maxima whenever ${\mathcal {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 aforementioned co-occurrence between transitivity peaks and modularity troughs, it is remarkable that this observation applies to all EP El Niño episodes, while there is a single El Niño event (1968/69) without transitivity peak (and hence classified at CP type according to [27] consistent with several previous studies referenced therein) but exhibiting a marked drop in modularity. This could point to either a specific type of ENSO event characterized by the emergence of teleconnectivity without markedly densified regional connectivity in the tropical Pacific Ocean, or the emergence of teleconnectivity unrelated to ENSO. We leave a more detailed analysis of this specific CP El Niño episode as a subject of future work.

4.1.2 Global topological network characteristics and ENSO transitional complexity

In addition to ENSO diversity, i.e. the distinction of El Niño and La Niña episodes into different flavors, the associated transitional complexity has recently attracted considerable interest among the scientific community [84,85,86]. The latter aspect refers not to individual ENSO events, but rather their successions, and gives rise to distinguishing “isolated” El Niño and La Niña episodes not followed by another event in the subsequent year, successions of El Niño and La Niña episodes with an abrupt reversal in tropical Pacific SST anomalies in between, and even multi-year events of the same type.

Unfortunately, the overall number of El Niño and La Niña events (of both flavors) contained in our data set does not allow for a systematic quantitative statistical investigation of the network characteristics associated with all possible variations of transitional behaviors. Indeed, because of the limited time coverage of reliable observational and reanalysis data sets, the transitional complexity of ENSO has been commonly studied based on long-term climate model simulations [86]. Further exploring this aspect by means of functional climate networks provides an interesting subject for future research, but is beyond the scope of the present work.

Fig. 3

Global maps showing composites of (a, c, e, g, i) degree $k_i$ and (b, d, f, h, j) average link distance $d_i$ for different types of ENSO phases: (a, b) EP El Niño, (c, d) CP El Niño, (e, f) EP La Niña, (g, h) CP La Niña and (i, j) all other periods. The corresponding classification of different years is summarized in Table 2. For constructing the composite maps, time windows corresponding to each type of event have been selected according to their midpoints coinciding with the end (31 December) of the calendar year given in Table 2

In the following, we will nevertheless take a closer look at the short-term temporal variability features of transitivity and modularity that arise in transitional complexity cases that can be identified in our reanalysis data set. Here, we will neglect transitions solely involving CP type events, since the latter have been shown above (as well as in [26, 27]) not to result in marked climate network transitivity (and in most cases also modularity) signatures. This qualitative analysis will help us to uncover (more subtle) differences related to the emergence and magnitude of the respective transitivity and modularity anomalies associated with EP type ENSO events exhibiting distinct transition sequences that go beyond the markedly negatively correlated interannual variability of both network measures.

1. Isolated EP El Niños (La Niñas) not followed by La Niña (El Niño) in the subsequent year. In our data set, we find three EP El Niños (1957/58, 1965/66 and 1982/83) and one EP La Niña (2008/09) that have not been followed by another ENSO event in the subsequent year. Commonly, such isolated EP events are characterized by large magnitudes of SST anomalies and a relatively large temporal persistence. In terms of our climate networks, all four events are characterized by two distinctive peaks in transitivity along with two simultaneous breakdowns in modularity, which peak at time windows whose midpoints approximately coincide with the beginning and end of the El Niño or La Niña conditions. In all four cases, the second transitivity peak and modularity trough are markedly weaker than the first one (see Fig. 2a and b).

Given the known seasonal profile of El Niño peaking around Christmas, it is remarkable that for all four isolated ENSO events, the ONI index has remained high (low) during a quite long period of time, resembling a single extended event, whereas the network measures rather seem to indicate two events—despite the absence of an immediately following opposite ENSO event in which case such a two-peak pattern could be expected (see below). To understand the apparently counter-intuitive and asymmetric behavior of ${\mathcal {T}}$ and ${\mathcal {Q}}$, we recall that ${\mathcal {T}}$ measures the presence of transitive structures, which most likely originate primarily from the highly synchronized regional dynamics in the ENSO region. Note again that longer links are commonly less likely to occur in a climate network than shorter ones, so that it is easier to generate transitive structures at smaller spatial scales than at larger ones. In contrast to this, ${\mathcal {Q}}$ captures changes in the hierarchy of connected patterns, which can be found across various possible spatial scales.

Specifically, maxima in ${\mathcal {T}}$ indicate the predominance of densely connected localized structures typical for strong EP El Niño and La Niña episodes as opposed to their CP counterparts [26], which is also consistently observed in Fig. 3a,e. In turn, the co-occurring minima in ${\mathcal {Q}}$ additionally highlight the relevance of the simultaneous emergence of teleconnectivity, i.e., long-range connections originating from the eastern tropical Pacific (Fig. 3b,f), turning the resulting climate networks into less modular structures by connecting otherwise mostly disconnected regions of the globe. Regarding the latter, the corresponding spatial patterns indicate that these teleconnections are mostly inner-tropical between the three relevant ocean basins, but partly also connect tropics and extratropics consistent with state of the art findings of climate research [41]. Further identifying and quantifying these teleconnections in our network representations would require more detailed studies on the link placement in the associated networks by properly regionalized measures, which is beyond the scope of the present work, but provides an interesting aspect to be further investigated in possible follow-up studies. Taken together, the above discussion underlines that both ${\mathcal {T}}$ and ${\mathcal {Q}}$ actually capture different aspects of network organization that provide complementary information.

The emergence of two distinct anomalies of ${\mathcal {T}}$ and ${\mathcal {Q}}$ along with “isolated” EP events can be understood when taking the temporal extent and evolution of those events into account. In the initial development phase of El Niño or La Niña, consistent SST anomalies emerge and subsequently increase over a vast part of the tropical Pacific Ocean. Given a sufficient temporal overlap between the time window for which the climate network is constructed and this development phase, the large-scale trends elevate the local correlations, along with the emergence of long-distance cross-correlations. The first peak (trough) of transitivity (modularity) can therefore be interpreted as a hallmark of this highly synchronous development phase. After the El Niño or La Niña conditions have been fully developed (in late boreal autumn), the SST anomalies and teleconnections remain relatively stable over a certain period of time. If this phase is sufficiently long, correlations between grid points in the corresponding time window used for network construction mainly reflect “background” fluctuations rather than common trends, leading to a reduction of the magnitude of those correlations and, hence, a tendency of transitivity and modularity to also return to their “background” values. Finally, the retreat of El Niño or La Niña conditions is again associated with common large-scale temperature trends (and, hence, elevated correlations), which may, however, appear less synchronously than during the development phase. As a consequence, the second transitivity maxima and modularity minima are commonly less well expressed than the first one coinciding with the development phase.

In any case, it is notable that the emergence of two directing transitivity and modularity anomalies along with isolated EP-type ENSO phases is fostered by the selected temporal window width of one year used throughout this work. While using shorter time windows might further magnify the corresponding effects, they also bring about substantial seasonal variability (not shown) that would undermine a proper interpretability of the observed time evolution of network properties.

2. EP El Niños (La Niñas) followed by La Niñas (El Niños) in the subsequent year. A markedly different situation with a however similar evolution of the two topological network characteristics ${\mathcal {T}}$ and ${\mathcal {Q}}$ arises in the context of immediate successions between El Niño and La Niña phases. In our data set, we find four successions of EP El Niños directly followed by EP La Niñas (1973–75, 1986–89, 1997–99 and 2009–11) and one of the opposite order (1975–77). All five cases are again characterized by two marked transitivity peaks associated with modularity minima, which approximately occur for those time windows with midpoints coinciding with the beginning of the El Niño or La Niña phase. While this is very similar to the case of isolated events, we find as a distinctive difference that for two-event successions, the second maximum/minimum commonly appears with similar or even larger magnitude than the first one. This can be expected because of the very rapid change between strongly positive and strongly negative tropical Pacific SST anomalies within about half a year (i.e. less than the window width used for network construction), implying very strong common trends along with the retreat phase of the first and following development phase of the second (opposite) ENSO event.

In the extreme case of the shift from (very strong) El Niño to La Niña conditions in summer 1998, a particularly fast reorganization of the global SAT field took place. The latter transition is reflected by some negative anomaly of ${\mathcal {T}}$ with respect to the baseline value for time windows centered in boreal autumn 1998, 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 ${\mathcal {Q}}$. This indicates that the surface air temperature field has been spatio-temporally re-organized primarily ín the ENSO region after the onset of La Niña, while teleconnections have already lost their relative importance as compared to the central periods of El Niño or La Niña phases as visible in terms of monotonically increasing modularity values in Fig. 2.

3. Multi-year El Niños (La Niñas). The last notable case of transitional complexity is the multi-year La Niña of 1954-56, which was characterized by persistent cold anomalies in the tropical Pacific Ocean over two boreal winter seasons. Such behavior can be seen as the extreme case of persistent “isolated” ENSO events or sequences of two ENSO events of the same type (e.g. the two El Niño years 1968/69 and 1969/70). Strictly speaking, the 1986–88 El Niño episode also covered two ENSO seasons (with the pronounced transitivity and modularity responses occurring during the first one), but was followed by an immediate La Niña period in 1988/89 and has therefore already been described above. It is notable that two-year events may be composed of a succession of EP and CP type events, accordingly showing only a single transitivity peak and modularity trough or even no significant anomaly. Since we do not have more events of this type on record, it is not possible to draw any further conclusions on the behavior in case of multi-year ENSO events from the present analysis.

Discussion and outlook. Our detailed analysis of the two global topological network characteristics transitivity and modularity suggests that not only the ENSO diversity, but also the transitional complexity of ENSO events contributes markedly to the emerging climate network features. While network modularity and transitivity generally evolve oppositely at interannual time scales, both provide complementary viewpoints on the emergence of elevated regional versus inter-regional connectivity. High transitivity commonly coincides with the temporary appearance of densely connected structures in the functional climate network constructed from global SAT anomalies, which are typically well localized in space [26]. 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. Regarding the differential effects of different ENSO episodes’ flavors and transitional behaviors, we have identified characteristic features in the emerging network structures that go beyond previously described differences due to ENSO diversity (i.e. EP versus CP events). Due to the very limited number of events contained in the considered reanalysis data set, fully exploring the reported findings and supplementing them by more quantitative statistical analysis will require larger data sets, most likely originating from long (control) runs or large ensembles of state of the art coupled climate model experiments. Corresponding analyses, which we outline here as important avenues for follow-up research, may address, among others, two important aspects not investigated here:

• 1. How characteristic and statistically significant are the qualitative findings described above on the different transitivity and modularity signatures along with ENSO diversity and transitional complexity, taking into account the underlying seasonality of the ENSO phenomenon controlling the timing of the corresponding network reorganization processes? Being provided with a sufficiently large sample size, this aspect could be addressed using composites and superposed epoch analysis, along with proper statistical significance testing.

• 2. How can we explain the fact that the largest ENSO magnitudes do not always coincide with the strongest anomalies in the climate network characteristics? In this context, it is important to note that the strength of an El Niño (or La Niña) event is described by the magnitude of warming/cooling (and, because of the seasonality of ENSO, also the associated rate of SST changes) in the tropical Pacific Ocean, but (potentially) also affects the spatial extent and persistence of those anomalies. By contrast, we do not necessarily have to expect a direct effect on the magnitude of changes in network properties resulting from spatial co-variability of those SST changes. The latter aspect should be rather affected by (i) the synchronicity of ENSO-related SAT changes across vast parts of the tropical Pacific, and (ii) by whether changes in the connectivity are more confined to the tropical Pacific or also trigger strong long-distance teleconnections. It appears reasonable that both can markedly depend on additional factors like the state of other climate variability modes like the Pacific Decadal Variability, Pacific Meridional Mode, Indian Ocean Dipole, Atlantic Nino, or Madden-Julian Oscillation, which operate on intraseasonal to interdecadal time scales and precondition the emergence of possible global-scale responses to El Niño and La Niña, but also the emergence of those ENSO episodes themselves. As a consequence, large El Niño and La Niña amplitudes may not be expected to coincide with the strongest changes in the network characteristics.

4.1.3 Global spatial network properties: average link distance

While the two above studied measures transitivity and modularity present key topological network characteristics, functional climate networks are systems embedded in geographical space. Thus, the spatial placement of nodes and links (which is not explicitly accounted for by topological characteristics) can play a pivotal role in network structure formation [26]. To address this aspect, we finally present the temporal evolution of the global average link distance $\left\langle \left\langle d\right\rangle \right\rangle$ 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 ${\mathcal {T}}$ and ${\mathcal {Q}}$. Yet, the general behavior of $\left\langle \left\langle d\right\rangle \right\rangle$ displays certain similarities with respect to that of network transitivity in the sense that ENSO-related peaks often, but not always co-occur in both measures, yet with marked time lags of up to several months. This lagged co-variability is clearly visible in the ENSO seasons 1957/58, 1965/66, 1973/74, 1982/83, 1997/98, 1998/99 and 2007/08 (with only 1982/83 showing two $\left\langle \left\langle d\right\rangle \right\rangle$ peaks along with an “isolated” El Niño event, while all other mentioned episodes rather display a single peak) and indicates that strong El Niño and La Niña episodes do not exclusively trigger short-range (localized) connectivity (high ${\mathcal {T}}$), but also global teleconnectivity (high $\left\langle \left\langle d\right\rangle \right\rangle$). This finding is in line with contemporary knowledge on the large-scale impacts of both types of ENSO phases [41, 87] and agrees well with previous qualitative results of [26] on the link distance distribution of global SAT networks.

It has to be noticed that considering a single value of the fixed link density of 0.5% may not be the optimal analysis strategy when facing a situation where the correlations within the ENSO regions and with other teleconnected regions are amplified by a different magnitude. In such cases, we may, for example, miss the effect of emerging teleconnectivity in the average link distance while being able to observe densified regional connectivity in the network transitivity. Accordingly, the results described here should always be interpreted as conditional with respect to the chosen link density. For the case of transitivity and modularity, Appendix A demonstrates that the overall interannual variability pattern is not affected qualitatively when increasing or decreasing the link density by a factor of 5. For the average link distance, which already presents “more noisy” interannual variability at the link density of 0.5% empirically determined as a reasonable choice in [26], the corresponding effect (not shown here for brevity) is however potentially more critical for the interpretability of the obtained results.

Beyond the aforementioned examples displaying a lagged co-occurrence between transitivity and average link distance peaks, we also find interesting cases where the link distance lacks any marked peak despite the presence of a strong EP type ENSO episode, particularly along with the ENSO seasons 1972/73, 1986–88 and 1988/89, indicating the absence of long-distance teleconnections (at least for the given correlation threshold). In turn, in various cases (1965/66, 1970/71, 1973/74, 1982/83, 1988/89 and 1997/98) we find the corresponding El Niño or La Niña periods being preceded by a drop in the average link distance. We interpret this finding such that in those cases, the emergence of strong regional (short-range) connectivity in the tropical Pacific actually preceded the formation of (long-distance) teleconnections, which implies first a drop in average link distances followed by a later peak. This could indicate that the emerging SST anomalies in the tropical Pacific trigger large-scale teleconnective effects to a greater extent than being themselves triggered by large-scale teleconnections. In fact, recent research especially on tropical basin interactions [88,89,90], but also on preconditioning of ENSO by North Pacific SST anomalies [91], points to the fact that the emergence of ENSO episodes in the tropical Pacific climate is causally coupled with key variability modes in other regions in both ways. In this context, the observed behavior of the average link distance of evolving SAT networks is well in line with those recent findings.

Unfortunately, the more “noisy” behavior of the average link distance as compared to network transitivity and modularity does not allow us making more conclusive statements at this point. One way to further pursue into this direction could involve systematically varying the link density to check if there is some optimum characterized by a larger “signal/noise ratio”, as well as differentiating more systematically between links connecting nodes within and between specific regions that could be related to previously reported ENSO teleconnections. A corresponding approach might also provide a way to further explain two additional observations associated with the variability of $\left\langle \left\langle d\right\rangle \right\rangle$ that have not yet been discussed above.

One one hand, we find an unprecedented isolated peak of $\left\langle \left\langle d\right\rangle \right\rangle$ in 1968 (indicating emerging teleconnectivity) unrelated to any obvious ENSO episode previously emerging (note that the anomalous modularity signature of the 1968/69 CP El Niño arises considerably later, not earlier, than the unexplained peak in the average link distance). This could indicate a case where teleconnections developed along with “sub-critical” ENSO conditions not leading to the formation of El Niño or La Niña, the specific reasons for which would deserve further more detailed studies.

On the other hand, two successive peaks of $\left\langle \left\langle d\right\rangle \right\rangle$ emerged in early 1963 and early 1964, i.e. before and after the short-lived 1963/64 El Niño episode, which coincides with an unprecedented transitivity peak and modularity trough found only after the retreat of the El Niño related tropical Pacific SST anomalies, whereas the expected transitivity peak and modularity trough at the beginning of this “isolated” event had been largely suppressed. The latter could indicate a reversal of the more common behavior along with strong EP El Niño phases as discussed earlier in this work, corresponding to the following “mechanistic” interpretation. First, starting from elevated teleconnectivity (large average link distance for time windows centered in early 1963), the emergence of El Niño conditions in the tropical Pacific Ocean may have been remotely forced instead of internally generated. One potential external triggering factor could have been the strong Mount Agung eruption starting in early 1963 [61]. Second, these remotely forced regional SST anomalies developed in a relatively heterogeneous/non-synchronous manner (suppressed transitivity and modularity) and did not lead to very strong and persistent El Niño conditions. Finally, the retreat of those “unstable” tropical Pacific SST anomalies occurred in a more coordinated manner (large transitivity/low modularity) than their emergence and was associated with the re-emergence of (delayed) long-distance connections in the climate network induced by the successive cooling of the tropical Pacific, which however did not reach the La Niña “threshold”. To this end, this possible interpretation necessarily remains speculative and requires further assessment.

From the results discussed above, we conclude that to distinguish globally influential ENSO events from episodes of minor (or more regional) relevance, a combination of transitivity with modularity and average link distance can be useful, taking a holistic perspective in studying the differential imprints of different types and successions of ENSO phases. We will recall this strategy when discussing the effects of volcanic eruptions on the network organization at a global scale.

4.1.4 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 patterns of the generated functional climate networks. Specifically, following recent observations that climate network properties distinguish between the EP and CP flavors of both El Niño and La Niña [27], 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 have been classified as showing either of these situations (see Table 2).

The left panels of Fig. 3 display the respective mean degree fields for the different types of ENSO periods. As expected, 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 behavior is expected from previous studies [27]. 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 flavors of El Niño and La Niña actually generate additional connections both within and out of the tropical Pacific that span relatively large distances, but possibly to different degrees and with different spatial teleconnectivity patterns [42, 87].

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 identified 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 [42] and, thus, via airflow in higher altitudes rather than near-surface atmospheric circulation. With a relative shift in the position of main tropical Pacific SST anomalies between EP and CP episodes, the regions exhibiting up-/downwelling airmasses are also shifted east-/westward. Hence, the anomalies associated with CP phases may affect regions more to the west than those observed along with EP phases, e.g. the Indian Ocean region rather than the maritime continent [87]). While accordingly elevated degree and average link distance values in the central Indian Ocean are hardly visible in Fig. 3, the corresponding results of Wiedermann et al. [27] who used a weighted network representation indicated slightly enhanced connectivity in that region along with CP episodes. In this context, it is notable that the choice of a link density of only 0.5% in our present work was designated to filter out only the strongest links (largest absolute correlation values) among different locations, thereby bearing some risk of missing elevated but still sub-threshold values in specific regions. Systematically varying the link density could provide a way to better highlight additional teleconnections not visible in Fig. 3, but exceeds the limitations of the present manuscript.

Moreover, it has to be noted that our present analysis is based on cross-correlations only. The values of this statistical similarity measure 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 [92, 93]. Accounting for this effect in terms of replacing the correlation values by associated significance levels in the network generation step [94] 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 [95].

The teleconnection with Western Africa spans a rather large spatial distance (about one third of the globe). In this context, Joly and Voldoire [96] noted that “a significant part of the West African monsoon (WAM) inter-annual 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 (which typically have larger absolute SST anomalies than their CP counterparts).

In general, climate variability within the tropics is typically more likely to exhibit strong correlations than between tropical and extratropical regions, 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 could merely be data and/or interpolation artifacts resulting from our remapping procedure rather than indications of actual teleconnections. Note that there exist polar-tropical teleconnections (and even cross-equatorial teleconnections) commonly mediated via the stratosphere [97, 98], which we however do not expect to be observable in lag-zero co-variability. Further examining the causes of elevated high-latitude connectivity requires better understanding the origins of the underlying long-distance links, which should be attempted in future work.

Related to the aspect of extratropical ENSO related teleconnections, we emphasize that the employed methodology of functional climate network analysis provides a new avenue for studying tropical basin interactions, i.e. processes linking coupled ocean-atmosphere variability in the three tropical ocean basins with each other via atmospheric (e.g. Walker circulation) and oceanic (e.g. regulation of cross-basin water transport through the Indonesian Throughflow). Recent work has provided multiple evidence for such linkages, and the elevated degree and average link distance over the equatorial Atlantic Ocean and maritime continent may be further assessed in the context of possible linkages between ENSO, parts of the Indian Ocean Dipole (IOD), and the Atlantic Niño and Atlantic Meridional Mode, respectively. Indeed, the Atlantic Niño has been discussed as a possible driver of El Niño [88] (including the use of a different type of complex network approach [99]), while the interrelationship between ENSO and IOD [89, 90] is still far from being conclusively understood. To this end, we suggest that a more detailed picture on those tropical basin interactions could be obtained using regionalized network properties (like being employed for the ENSO region only in the following) and/or applying the concept of coupled climate networks [23], along with replacing instantaneous by lagged correlations to also account for the common time lag between relevant changes in the atmosphere-ocean coupled variability in the three tropical basins. Further analysis on this aspect is indeed planned as an integral part of our future work.

4.1.5 Regionalized network characteristics

Fig. 4

Time series of different regionalized climate network properties introduced in Sect. 3.3.1. Background colours and time axis are the same as in Fig. 2

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 regional structure, we would expect a strong increase of $k_i$ but only a 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. In this regard, the information provided by the regionalized degree and average link distance for the three considered regions is partially redundant as well. 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 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. Hence, it might be more relevant to look at differences between the respective regionalized degree and average link distance variations. Another valid concern in this context is the strong inter-connectivity between the different considered regions, which may dominate especially the regionalized degree in comparison with extra-Pacific or tropical-extratropical teleconnections. However, the latter would be partially indicated by our approach, so that studying both types of network characteristics together may still provide valuable information. A methodological alternative would be decomposing the entire globe into different regions of interest (e.g., different parts of the tropics and extratropics) and use the framework of coupled network analysis [23] for studying the presence or absence of linkages between these regions in a more detailed way. We outline further investigations in this direction as a subject for future work.

Examining the time evolution of all six regionalized network measures in some detail, it is notable that between 1978 and 1982 (coinciding with the beginning of satellite measurements assimilated in the reanalysis data set, i.e., a source of heterogeneity in the data), 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 spatio-temporal 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.

4.2 Volcanic eruptions

Besides ENSO variability, strong volcanic eruptions have been identified as additional causes of marked disruptions in climate network properties in earlier studies [26]. 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 $\left\langle \left\langle d\right\rangle \right\rangle$, we find a marked difference. Specifically, for ENSO-related climate disruptions, peaks in ${\mathcal {T}}$ have been associated with both simultaneous drops and delayed peaks in $\left\langle \left\langle d\right\rangle \right\rangle$, indicating the presence of strong localized connectivity structures in the tropical Pacific and the excitation of teleconnections, respectively. By contrast, in 1992/93 we find a persistent positive transitivity anomaly associated with an extended through in modularity, but only a persistent drop (without later peak) in the average link distance. Notably, the corresponding period corresponds to the time windows where the tropospheric cooling effects following the Mount Pinatubo eruption should have taken their maximum [64]. Hence, unlike for ENSO-related disruptions, the peak in transitivity together with a simultaneous drop in $\left\langle \left\langle d\right\rangle \right\rangle$ indicate that the effects of the volcanic eruptions have rather been regionally confined (i.e., we observe the emergence of strong regional connectivity due to large-scale coherent trends, but no excitation of relevant long-distance correlations). The latter is consistent with the presence 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 distance in combination 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).

Fig. 5

Time series of regional mean degree and associated degree field (excluding the ENSO-big region indicated by white color) for the three strongest 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^{\circ }$ 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 colors indicating the corresponding ENSO strength are the same as in Fig. 2

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 [55]. Along with 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 is 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 large-scale cooling effect (also markedly affecting the global mean temperature) 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 the regionalized degree before the actual eruption, pointing to a possible interference with normal natural variability (specifically, approaching ENSO events that might have been (partially) triggered by the eruption [61]).

Unlike for 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, since El Chichon is located relatively close to the extended ENSO region, it is possible 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 [61]. In general, previous studies indicate that the El Chichon eruption caused a much weaker summer cooling than the Mount Agung eruption [100], which could also explain its absent signature in our analysis.

Taken together, our results regarding the climate network reorganization along with the largest three volcanic eruptions since the mid-20th century clearly indicate that strong enough eruptions can result in the emergence of densely connected network structures without additional long-ranging spatial correlations as typical for strong ENSO events. Thereby, while the previously considered network measure transitivity would mostly indicate the former effect (i.e., the emergence of spatially localized dense connectivity structures), combining this information with spatial characteristics as captured by the global average link distance $\left\langle \left\langle d\right\rangle \right\rangle$ allows us to discriminate large-scale climate disruptions regarding whether or not they have excited any relevant teleconnectivity. In this context, it is worth to recall again the differential climate effects of volcanic eruptions depending on the strength, location and season of the event. While the proposed network-based approach is potentially useful for categorizing and quantifying these effects, a specific process-based attribution would require consideration of additional information, i.e., the use of multiple relevant climatological fields instead of SAT anomalies alone as studied in this work. We outline corresponding follow-up studies as a prospective task for future research.

5 Conclusions and outlook

We have used functional climate networks constructed from spatial correlations of daily global surface air temperature (SAT) anomalies to further analyze the well-established emergence of distinctive spatio-temporal co-variability patterns associated with past El Niño and La Niña events as well as strong volcanic eruptions. Here, the main motivation was to obtain a better understanding of the climatological interpretation of different network characteristics, which is of key importance for guiding the selection of network properties to be used in future studies on more specific problems. In this context, by making use of the global network property of modularity (${\mathcal {Q}}$), we have found that the East Pacific flavors of El Niño and La Niña events lead to a global reconfiguration of SAT variations in terms of emerging teleconnections, whereas other climate disruptions have rather been characterized by regional changes only. The corresponding perspective is markedly distinct from that of previous studies of climate networks, which have mostly focused on network transitivity (${\mathcal {T}}$) as a discriminative statistics of different types of ENSO phase [26, 27]. The results obtained in the present work clearly indicate that ${\mathcal {T}}$—unlike ${\mathcal {Q}}$—may be primarily sensitive to the presence of common regional trends (in case of ENSO, mostly in the eastern tropical Pacific) and less indicative of the emergence of actual teleconnections spanning vast parts of the globe. By 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 co-variability 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 [95, 96]. 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 SAT globally, its dominating effect was rather regional (i.e., it did not trigger marked additional 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. [27] presenting 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.

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 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, leading to generally weaker climatic effects than for tropical eruptions. 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 [101].

The fact that ENSO does not only have a strong regional effect on SAT in terms of coherent SAT trends in the tropical Pacific associated with a spatially confined increase of network connectivity [26], but also dynamically reconfigures climate variability globally by triggering teleconnections especially with other tropical regions, has been previously found in other studies utilizing more classical analysis techniques. In this regard, one possible mechanism could involve the modulation of monsoons by strong El Niño and/or La Niña periods, which itself could be further affected by volcanic eruptions [63, 102]. 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.

To this end, we are confident that the combination of different structural and spatial climate network properties provides a prospective tool for making the qualitative statements discussed in this work more quantitative and, thus, discriminating further between large-scale climate effects due to regional trends and such originating from modifications of teleconnectivity. Notably, network properties can serve as a unique toolbox for disentangling both commonly closely interconnected aspects, and hence offer some interesting nonlinear research methodology for further research focusing on specific climatological processes, e.g., in combination with ensembles of simulations of state-of-the-art general circulation models. By means of systematic applications of these tools on such vast amounts of high-quality (and at least partially “controllable”) climate data, it is likely that important new insights may be gained that (as in case of other successful applications of complex network theory in climatology as mentioned in the Introduction) classical statistical analysis techniques might not be able to provide.

Fig. 6

Network transitivity (a) and modularity (b) evolution for a link density of 0.1%

Fig. 7

Same as Fig. 6 for a link density of 2.5%

Fig. 8

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 WalkTrap 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

Appendices

Appendix A: Dependence of transitivity and modularity on the link density

The selection of a link density of 0.5% in the main body of this paper has been based on detailed investigations in previous studies [26, 27]. In general, it is a heuristic choice that has been found to provide a reasonable trade-off in the context of making the network characteristics relatively sensitive to the emergence of “localized structures” in the climate network. This does not mean that other choices are generally worse, but for the purpose of this study (tracing subtle changes in the global spatio-temporal correlation pattern of surface air temperatures), too high link densities might blur the effects of interest while too low densities could make the network decompose into disconnected components, thereby prohibiting the application of the selected network characteristics.

To illustrate the choice of 0.5%, as well as the qualitative robustness of some of our main findings, we have repeated our analysis with two different link densities of 0.1% and 2.5%, respectively, and did not find strong differences to our previous analysis (see Figs. 6 and 7 for some results).

Appendix B: Comparison of different modularity estimation algorithms

In Fig. 8, 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 [111], infomap [112], label propagation [113], leading eigenvector [114] and WalkTrap [76]. Further details motivating the choice of WalkTrap as a reference algorithm in the body of this paper are provided in the figure caption.