Articles | Volume 26, issue 3
Nonlin. Processes Geophys., 26, 251–266, 2019

Special issue: Centennial issue on nonlinear geophysics: accomplishments...

Nonlin. Processes Geophys., 26, 251–266, 2019

Research article 15 Aug 2019

Research article | 15 Aug 2019

Unravelling the spatial diversity of Indian precipitation teleconnections via a non-linear multi-scale approach

Unravelling the spatial diversity of Indian precipitation teleconnections via a non-linear multi-scale approach
Jürgen Kurths1,2,3,,**, Ankit Agarwal1,2,4,, Roopam Shukla1, Norbert Marwan1, Maheswaran Rathinasamy5, Levke Caesar1,6, Raghavan Krishnan7, and Bruno Merz2,4 Jürgen Kurths et al.
  • 1Potsdam Institute for Climate Impact Research (PIK), Member of the Leibniz Association, Telegrafenberg, Potsdam, Germany
  • 2Institute of Earth and Environmental Science, University of Potsdam, Potsdam, Germany
  • 3Institute of Physics, Humboldt Universität zu Berlin, Germany
  • 4GFZ German Research Centre for Geosciences, Section 5.4: Hydrology, Telegrafenberg, Potsdam, Germany
  • 5Dept. of Civil Engineering, MVGR College of Engineering, Vizianagaram, India
  • 6Institute of Physics and Astronomy, University of Potsdam, Potsdam, Germany
  • 7Indian Institute of Tropical Meteorology, Pune, India
  • These authors contributed equally to this work.
  • **Invited contribution by Jürgen Kurths, recipient of the EGU Lewis Fry Richardson Medal 2013.

Correspondence: Ankit Agarwal (,,


A better understanding of precipitation dynamics in the Indian subcontinent is required since India's society depends heavily on reliable monsoon forecasts. We introduce a non-linear, multiscale approach, based on wavelets and event synchronization, for unravelling teleconnection influences on precipitation. We consider those climate patterns with the highest relevance for Indian precipitation. Our results suggest significant influences which are not well captured by only the wavelet coherence analysis, the state-of-the-art method in understanding linkages at multiple timescales. We find substantial variation across India and across timescales. In particular, El Niño–Southern Oscillation (ENSO) and the Indian Ocean Dipole (IOD) mainly influence precipitation in the south-east at interannual and decadal scales, respectively, whereas the North Atlantic Oscillation (NAO) has a strong connection to precipitation, particularly in the northern regions. The effect of the Pacific Decadal Oscillation (PDO) stretches across the whole country, whereas the Atlantic Multidecadal Oscillation (AMO) influences precipitation particularly in the central arid and semi-arid regions. The proposed method provides a powerful approach for capturing the dynamics of precipitation and, hence, helps improve precipitation forecasting.

1 Introduction

Understanding the spatial patterns, frequency and intensity of precipitation in the Indian subcontinent is an active area of research due to its essential impact on life and property (Pai et al., 2015). The Indian monsoon is the pulse and lifeline of over 1 billion people, and the socio-economic development in this part of the world heavily depends on reliable predictions of the monsoon (Goswami and Krishnan, 2013; Shukla et al., 2018).

Numerous studies have emphasized the importance of understanding the influence of large-scale climatic patterns on precipitation for improving forecast accuracy (Feng et al., 2016), and therefore many studies have analysed the relationship between precipitation and climatic patterns for India. This research has shown that the relevant patterns are the El Niño–Southern Oscillation (ENSO) (Kumar et al., 2006; Mokhov et al., 2012), the Indian Ocean Dipole (IOD) (Behera et al., 1999; Krishnan and Swapna, 2009), the North Atlantic Oscillation (NAO) (Bharath and Srinivas, 2015; Feliks et al., 2013), the Pacific Decadal Oscillation (PDO) (Dong, 2016; Krishnan and Sugi, 2003), and the Atlantic Multidecadal Oscillation (AMO) (Goswami et al., 2006; Krishnamurthy and Krishnamurthy, 2016).

Over the years, linkages between climatic patterns and precipitation have been investigated by a range of statistical methods, such as correlation (Abid et al., 2018), principal component analysis (Luterbacher et al., 2006), empirical orthogonal functions (Hannachi et al., 2007), and regression and canonical analysis (Xoplaki et al., 2004). However, all these methods are limited in capturing the scale-specific feedbacks and interactions between the long-range climatic patterns and precipitation. Such information is very crucial since in climatic systems energy is stored and transported differently on different temporal scales, resulting from interactions of intertwined sub-components across a wide range of scales (Miralles et al., 2014; Peters et al., 2007). Multiscale interactions have therefore received extensive attention in the field of climate dynamics (Peters et al., 2007; Steinhaeuser et al., 2012) and have been proposed as a mechanism for triggering extreme events (Agarwal et al., 2018b; Okin et al., 2009; Paluš, 2014; Peters et al., 2004) and abrupt transitions (Peters et al., 2007). This holds the promise of better understanding the system dynamics compared to analysing processes at one timescale only.

In recent decades, wavelet coherence has become the state-of-the-art method for studying the influence of climatic patterns on precipitation at different temporal scales. For example, Ouachani et al. (2013) investigated the multiscale linear relationship between the Mediterranean region (northern Africa) and large-scale climatic patterns such as ENSO, NAO and PDO. The study reported a strong correlation between ENSO and precipitation series at a lag of 2 years. The study further reported that the influence of ENSO on precipitation was stronger compared to other climatic modes considered in this particular study. Coherently, Tan et al. (2016) analysed the relations between Canadian precipitation and different global climate indices. Similar studies using wavelet coherence also reported in other parts of the world (Agarwal et al., 2016; Araghi et al., 2017; Hu and Si, 2016; Tan et al., 2016), though all such studies based on wavelet coherence were illuminating and contributed significantly to our existing knowledge of the climate. However, there remains a large scope for advancement, in particular in capturing the non-linear scale-specific interactions between climate patterns and Indian precipitation.

To capture such non-linear scale-specific interactions, recently event synchronization (ES) has emerged as a powerful similarity measure (Agarwal et al., 2019a; Mitra et al., 2017; Ozturk et al., 2018; Quiroga et al., 2002) because ES automatically classifies pairs of events arising at two locations as temporally close (and, thus, possibly statistically – or even dynamically – interrelated) without the necessity of selecting an additional parameter in terms of a fixed tolerable delay between these events. Also, ES is a robust measure to study the interrelationship between a series of non-Gaussian data or data with heavy tails (Agarwal, 2019). These intrinsic features of ES are advantageous in climate in general and for quantifying interactions between climatic patterns and precipitation in particular since the time delay between such patterns (for e.g. ENSO) and their effect on precipitation is tedious to quantify beforehand.

We, therefore, decided to use ES to quantify the (possibly non-linear) linkages between large-scale climatic patterns and precipitation across India. More specifically, we analyse the linkages between the 95th percentile extreme events, extracted from gridded Indian precipitation data at monthly resolution, and the climate patterns ENSO, IOD, NAO, PDO, and AMO which have been shown to be of significant relevance for precipitation in India. We combine ES with the wavelet transform, as proposed recently (Agarwal et al., 2017). This combination, termed MSES (multiscale event synchronization), allows non-linear connections between time series to be studied at different temporal scales. To consider the spatial variation across India, we sub-divide India into homogeneous regions that share rather similar precipitation characteristics and identify a representative grid cell for each region. The homogenous regions and the representative grid cells are obtained using the concept of a complex networks approach (Agarwal et al., 2018a), and the resultant network is referred to as a climate network (Agarwal et al., 2019a; Boers et al., 2019; Ekhtiari et al., 2019; Tsonis et al., 2006).

The novelty of this study is the integration of (1) the non-linear method for quantifying the linkages between large-scale climate patterns and the climate network (precipitation) in India at (2) multiple timescales, considering (3) the spatial variation of these linkages. To our knowledge, this combination (non-linear–multiple timescales–spatial variation) has not yet been implemented, neither for India nor for any other region. We argue that it allows the spatio-temporal diversity of Indian precipitation teleconnections to be unravelled, offering a compelling perspective for capturing the dynamics of precipitation and improving precipitation forecasting.

2 Study area and data

2.1 Study area

Our study area is the Indian subcontinent, which shows a significant variation in climate characteristics. India extends over an area of 3 287 263 km2. Its climate regimes are classified as arid (north-western India), semi-arid (northern lowlands and central peninsular India), humid (coastal lowlands, south-western and north-eastern highlands) and alpine (Himalayan mountains in the north). The spatio-temporal variation of precipitation, as well as temperature, is significant over the country (Bharath and Srinivas, 2015). The entire country receives 80 % of its total precipitation during the south-western monsoon, from June to September (Bharath and Srinivas, 2015). During the north-eastern monsoon (October to December), the precipitation is considerable but is confined to the south-eastern part of the country.

2.2 Gridded precipitation data

We use the high-resolution (0.25×0.25) monthly gridded precipitation dataset for the period 1951–2013, developed by the Indian Meteorological Department (IMD) for the spatial domain of 66.5 to 100 E and 6.5 to 38.5 N, covering the mainland region of India (Pai et al., 2014). The gridded data have been generated from the observations of 6995 gauging stations across India (Pai et al., 2014). The dataset captures well the spatial distribution of precipitation over the country. For our study, out of 17 415 grid cells, 4631 cells lying inside the boundaries of India were identified.

2.3 Time series of global and regional climate indices

To understand the linkages between climate patterns and precipitation, we use time series of global and regional climate indices for the same period, i.e. 1951–2013. We have selected those indices for which earlier studies have shown a relation to Indian precipitation. The selected climate indices and the respective studies are ENSO (Mokhov et al., 2012), IOD (Ashok et al., 2001), NAO (Bharath and Srinivas, 2015; Feliks et al., 2013), PDO (Dong, 2016; Krishnan and Sugi, 2003) and AMO (Goswami et al., 2006; Krishnamurthy and Krishnamurthy, 2016). For detailed information on these climate indices and the data sources, we refer the reader to (last access: 26 May 2018).

3 Methodology

To investigate the non-linear, multiscale linkages between climate patterns and precipitation, we propose an analysis based on combining network reconstruction, community detection, wavelet transformation, and event synchronization (Fig. 1). First, we construct a precipitation network of the precipitation dataset using event synchronization. We further pool grid cells with similar precipitation characteristics into homogenous regions and identify a representative grid cell for each region as proposed by Agarwal et al. (2018a). The linkages between the precipitation time series of the representative cells and the teleconnection indices are analysed by the MSES method developed by Agarwal et al. (2017). Finally, the proposed methodology is compared to the state-of-the-art wavelet coherence analysis (WCA).

Figure 1Schematic of the methodology to investigate the linkages between climate patterns and precipitation. (“M” stands for method and “R” for result.)

3.1 Event synchronization and network construction

ES measures the non-linear synchronization of point processes (Quiroga et al., 2002). Each grid cell of the precipitation dataset serves as a network node, while the precipitation estimates at each cell provide the time series for that node. Following Agarwal et al. (2018a), we define heavy precipitation events at each node as events with precipitation larger than the 95th percentile at that grid cell. The 95th percentile threshold for event selection is globally accepted as a tradeoff between a sufficient number of events and a high threshold value. Then ES is used to calculate the strength of synchronization (Q) between all possible pairs of grid cells. ES has advantages over other time-delayed correlation techniques (e.g. Pearson lag correlation), as it uses a dynamic (not fixed) time delay (Agarwal et al., 2017). The latter refers to a time delay that is adjusted according to the two time series being compared, which allows its application to different situations. Another advantage of ES is that it can be applied to non-Gaussian data. Having its roots in neuroscience, ES only considers events beyond a threshold and ignores the absolute magnitude of events, which could be a challenge to incorporate in future work.

Here, we define an event when a value in the signal x(t) (or y(t)) exceeds a threshold (selected by a α percentile). Events in x(t) and y(t) occurring at time tlx and tmy, where l=1,2,3,4Sx, m=1,2,3,4Sy, are considered to be synchronized when they occur within a time lag ±τlmxy which is defined as in Agarwal et al. (2017).

(1) τ l m x y = min t l + 1 x - t l x , t l x - t l - 1 x , t m + 1 y - t m y , t m y - t m - 1 y } / 2

Sx and Sy are the total number of events (greater than threshold α) that occurred in the signals x(t) and y(t), respectively. This definition of the time lag helps to separate independent events. Then we count the number of times an event occurs in the signal x(t) after the maximum time lag τlmxy of an event that appears in the signal y(t) and vice versa, resulting in the quantities C(x|y) and C(y|x):

(2) C x | y = l = 1 S x m = 1 S y J x y and C y | x = l = 1 S x m = 1 S y J y x ,


(3) J x y = 1 if 0 < t l x - t m y < τ l m x y , 1 2 if t l x = t m y , 0 else .

From these quantities, we define a measure of the strength of event synchronization (Qxy) between x(t) and y(t) by

(4) Q x y = C x | y + C y | x ( S x - 2 ) ( S y - 2 ) .

Qxy is normalized to 0Qxy1. Qxy=1 refers to perfect synchronization between the signals x(t) and y(t). ES has been specifically designed to identify non-linear associations among event time series with varying lags between them.

A link between two grid cells is set up if their heavy precipitation occurrences are strongly synchronized, which we define as having a Q value greater than a predefined threshold (θx,yQ). A number of criteria have been proposed to generate an adjacency matrix from a similarity matrix, such as a fixed amount of link density (Agarwal, 2019) or fixed thresholds (Donges et al., 2009). Here, we consider a 5 % link density since it is a well-accepted criterion in general for the network construction. Also, the 95th percentile is a good trade-off between having a sufficient number of connections and capturing high synchronized connections.

We repeat the procedure for all possible pairs of nodes to construct a precipitation network with the adjacency matrix

(5) A x , y = 1 , if Q x , y θ x , y Q , 0 , else .

Here, θx,yQ=95th percentile is a chosen threshold, Ax,y=1 denotes a link between the xth and yth nodes and 0 denotes otherwise. The adjacency matrix represents the connections in the rainfall network. In this study, we use an undirected network, meaning we do not consider which of the two synchronized events happened first, in order to avoid the possibility of misleading directionalities of event occurrences between nodes that are topographically close to one another.

3.2 Community detection and the ZP approach

The linkages between climate indices and precipitation are evaluated on a regional scale. India is subdivided into homogeneous regions with similar characteristics of heavy precipitation events using the concept of complex networks (Agarwal et al., 2018). Several studies such as Agarwal et al. (2019b), Halverson and Fleming (2015), Lancichinetti and Fortunato (2009), Newman (2006), Sivakumar et al. (2015), and Tsonis et al. (2011) have reported superior performance of complex networks in identifying homogeneous regions compared to more traditional methods, such as the hierarchical clustering algorithm or the information-theoretic algorithm (Harenberg et al., 2014).

There exist several community detection approaches aiming at stratifying the nodes into communities in an optimal way (see Fortunato, 2010, for an extensive review). The question of which community detection algorithm should be used is difficult to answer. However, it has been found that the choice of the community detection algorithm has a small impact on the resultant communities in geophysical data science studies (Halverson and Fleming, 2015). In this study, we use the Louvain method which maximizes the modularity to find the optimal community structure in the network. The optimal community structure is a subdivision of the network into non-overlapping groups of nodes, which maximizes the number of within-group edges and minimizes the number of between-group edges (Blondel et al., 2008; Rubinov and Sporns, 2011).

Modularity is defined, besides a multiplicative constant, as the number of edges falling within groups minus the expected number in an equivalent network with edges placed at random. Positive modularity values suggest the presence of communities. Thus, one can search for community structures by looking for the network divisions that have positive, and preferably large, modularity values (Newman, 2004). Modularity (M) is calculated as

(6) M = 1 2 m x , y A x y - k x k y 2 m δ ( C x C y ) ,

where Axy represents the number of edges between x and y,kx=yAxy is the sum of the number of the edges (degree) attached to vertex x, and Cx is the community to which vertex x is assigned, and the δ function δ(u,v) is 1 if u=v and 0; otherwise, m=12xyAxy.

Equation (6) is solved using the two-step iterative algorithm proposed by Blondel et al. (2008), also known as the Louvain method. The first step consists in optimizing the modularity by permitting only a local modification of communities; in the second step, the communities identified are pooled to assemble a new network of communities. High-modularity networks are densely linked within communities but sparsely linked between communities. The algorithm stops when the highest modularity is achieved.

Further, for each community, we identify a representative grid cell using the ZP space approach, where Z is the within-module degree or Z score and P is the participation coefficient (Agarwal et al., 2018a). The within-module degree (Zx or Z score) is a within-community version of degree centrality (total number of links of any node) and shows how well a node is connected to other nodes in the same community. It is estimated as in Guimer and Amaral (2005).

(7) Z x = K x - K s x σ k s x ,

where Kx is the total number of links (degrees) of node x in the community sx, Ksx is the average degree of all nodes in the community sx, and σksx is the standard deviation of K in sx. Since two nodes with the same Z score may play different roles within the community, this measure is often combined with the participation coefficient Px.

The participation coefficient (Px) compares the number of links of node x to nodes in all communities with the number of links within its own community. We define the Px of node x as in Guimer and Amaral (2005).

(8) P x = 1 - s y = 1 N M k x s y k x 2 ,

where kxsy is the number of links of node x to nodes in community sy and kx is the total number of links (degrees) of node xNM representing the number of communities in the network. The participation coefficient of a node is therefore close to one if its links are uniformly distributed among all the communities and zero if its entire links are within its own community because in the latter case Kxsy=Kx and hence Px=0.

The cell with the highest number of intracommunity links is considered representative (Halverson and Fleming, 2015), based on the argument that this cell shows the strongest synchronization within the community. We expect its climatological properties, such as the linkage to large-scale climate patterns, to have the highest similarity to the properties of the other cells in the community. We could also use a composite, e.g. by normalizing the grid cell time series and defining the time series of the mean of the normalized series as representative. However, this definition would reduce the variability and could mask existing connections to climatic patterns.

3.3 Multiscale event synchronization

In this study, we use the MSES measure (Agarwal et al., 2017) that combines the wavelet transform and event synchronization to quantify the relationship between precipitation and climate indices. The multiscale event synchronization measure is based on a combination of wavelet transform and event synchronization (Sect. 3.1). The following subsections discuss briefly the wavelet transform and finally the methodology for MSES.

3.3.1 Wavelet analysis

Synthetically the temporal data series of any continuous geophysical variable is the superposition of variations occurring at different scales. Different physical processes drive these patterns, and a partitioning of the variability at different scales can help to isolate and characterize underlying processes (Agarwal et al., 2018b; Agarwal et al., 2019a). Wavelets have been successfully used to characterize the timescale of interactions across fluxes and physical drivers (Katul et al., 2001; Ding et al., 2013).

The wavelet transform of a signal decomposes it into a set of components with predefined central frequencies and spectral bandwidths. Here we use the maximal overlap discrete wavelet transform (MODWT) (Percival and Walden, 2000) because the orthogonal discrete wavelet transform (DWT) results in a pyramid of wavelet coefficients which does not contain the time synchronization of the events. Further, our experience with DWT suggests that the latter approach suffers from “shift sensitivity”, also known as “shift variance”, and is undesirable because it implies that DWT coefficients fail to distinguish between input-signal shifts (Maheswaran and Khosa, 2012). Even though the MODWT has large redundancy, it is shift-invariant, and this property renders the MODWT more suited to time series analysis.

MODWT decomposes the time series into different timescales or frequency components. The wavelet decomposition is realized using the two basis functions known as father wavelets and mother wavelets. Any function f(t) can be expressed in these basis functions and their scaled and translated versions as given in Eq. (9).

(9) f t = k s J , k φ J , k t + k d J , k ψ J , k t + k d J - 1 , k ψ J - 1 , k t + k d 1 , k ψ 1 , k t ,

where J is the number of multiresolution components (scales) and k is in the range of 1 to the number of the coefficient in the specified component. The coefficients sJ,k are the approximation coefficients and dJ,k,…, d1,k are the wavelet transform coefficients, while the functions φJ,k(t) and ψj,kt)|j=1,J,-1,J are the approximating wavelet function and detailed wavelet functions, respectively.

These basis functions are defined in terms of father and mother wavelets as follows:




where the scaling coefficients sJ,k capture the smooth trend of the time series at the coarse scale 2J, which are also called smooth coefficients; and the wavelet coefficients dj,k, also known as detail coefficients can detect deviations from the coarsest scale to the finest scale.

The original series f(t) can be reconstructed by the summing the detailed components and the smooth components.

(14) f t = S J , k + D J , k + D J - 1 , k + D 1 , k ,


(15) S J , k = k s J , k φ J , k t , D J , k = k d J , k ψ J , k t D 1 , k = k d 1 , k ψ 1 , k t .

Equation (14) defines a multiresolution analysis (MRA) of f(t); i.e. we express the series f(t) as the sum of a constant vector SJ and J other vectors Dj, j=1,,J, each of which contain a time series related to variations in f(t) at a certain scale. We refer to Dj as the jth-level wavelet detail. Figure 2 shows the MODWT decomposition of a sample signal up to seven scales resulting in seven detailed components (D1–D7) and one approximate component (S7).

Figure 2Scheme of multi-scale decomposition of signals using maximum overlap discrete wavelet transformation (MODWT). The relationship between signal Yt (blue), detailed component Dj (black), and approximate component Sj (red) is shown.


Let Yt represent a time series history of a geophysical process. In order to partition the variability of the process at different scales j= 1…J the signal Yt is transformed into the wavelet space which provides the required information at different scales. This is obtained by convolving Yt with a set of low-pass (g) and high-pass (h) filters. For instance, at each scale j, the MODWT applies a high-pass wavelet filter hj,l and a lower-pass filter gj,l of length (l) to the time series Y to, respectively, yield the wavelet coefficients Wj,t and Vj,t for every point t in the time series (Percival and Walden, 2000).

(16) W j , t = l = 0 L j - 1 h j Y t -lmod2N V j , t = l = 0 L j - 1 g j , l Y t -lmod2N

The Wj,t wavelet coefficients distinguish fluctuations in the time series of scale 2j−1, while the Vj,t coefficients provide information about the variations at scale 2j and higher. Let the maximum level of decomposition be j=J. This would result in a total J+1 series of wavelet coefficients, with Wj,t, j=1,2,3J, and one series of VJ,t.

Let us now define Dj which represents the time domain reconstruction of Wj. It represents the portion of Y attributable to scale j. Let SJ represent the time domain reconstruction of VJ. For the maximum level of decomposition, VJ has all of its elements equal to the sample mean of Y.

Therefore, we can write

(17) Y = j = 1 J D j + S J .

3.3.2 Stepwise procedure to estimate MSES

The MSES values between precipitation and climate indices are estimated in the following manner.

  • a.

    The climate indices and precipitation values at monthly resolution are decomposed into its various scale-specific components as proxies of the corresponding signal using the maximum overlap discrete wavelet transformation (MODWT). These components represent the features of the signal at different timescales. We limit the analysis to scale 7, i.e. 16 years, due to the distortion created by the boundary effects of the wavelet decomposition (Percival, 2008).

  • b.

    After fixing a 95 % threshold for each of the decomposed components of precipitation and climate indices, the event synchronization values are estimated. The 95 % threshold values are estimated for each scale component separately, ensuring a reliable estimation of the synchronization between the events.

  • c.

    The estimated ES values are considered significant if they are higher than the ones obtained from a significance test (Agarwal et al., 2017).

  • d.

    These steps (a–c) are repeated for all combinations of climate indices and precipitation for the different regions.

3.4 Significance test for similarity measure

To evaluate the statistical significance of the ES values, a surrogate test is used as proposed by Agarwal et al. (2017). We randomly reshuffle each time series 100 times (arbitrary number) but keep the distribution the same. The reshuffling will ensure that any potential synchronization between the even series will be destroyed and that they will be equivalent to independent random series. Then, for each pair of time series (rainfall and climate time series), we calculate the MSES values for the different scales. At each scale, the empirical test distribution of the 100 MSES values for the reshuffled time series is compared to the MSES values of the original time series. Using a 1 % significance level, we assume that synchronization cannot be explained by chance if the MSES value at a certain scale of the original time series is larger than the 99th percentile of the test distribution.

3.5 Testing of MSES on a synthetic dataset

In a previous study, we have tested the MSES measure with different synthetic time series and have shown the efficacy of the method (Agarwal et al., 2017). In this paper, we further test the method with synthetic datasets and compare the results with those of the traditional methods such as correlation analysis and wavelet coherence.

Consider two time series X and Y (of length 1000) as defined by


where e(t) and H(t) denote a white noise process N(0,1) and a random binary series with values 1 or 0. H(t) represents the aperiodic extreme events in the given time series.

Figure 3 shows the plot of X and Y with respect to time. It can be observed that X and Y are time series with two distinct frequencies but have a lagged relationship induced by H(t). In this example, we have considered a=4. The zero lag correlation coefficient between X and Y can be estimated as −0.02 and the Pearson lag correlation is found to be 0.3, both showing no significant correlation. However, as expected, the MSES given by Agarwal et al. (2017) was estimated to be 0.9375 at scale 1, revealing the underlying synchronization between the two series in this scale only.

Figure 3Test signals X(t) (a) and Y(t) (b) with two distinct frequencies having a lagged relationship induced by H(t) used to explain the methodology.


3.6 Wavelet coherence analysis (WCA)

We benchmark the MSES results against the WCA because wavelet coherence is the state-of-the-art method in evaluating linkages between hydroclimatological variables at multiple timescales (Peters et al., 2004; Tan et al., 2016). We use the Grinsted Toolbox (Grinsted et al., 2004) for calculating the WCA between precipitation of the representative grid cells and the climatic indices. The wavelet coherence between time series X{Xt} and {Yt} was defined by Torrence and Compo (1998) as

(20) R 2 j , t = | ς ( j - 1 W x y ( j , t ) | ς ( j - 1 | W x j , t | 2 ) ς ( j - 1 W y j , t 2 ) .

ζ is a smoothing operator and can be written as ςW=ςscaleςtimeWj,t. Wxy represents the cross-wavelet coefficient between X and Y. Wx (j,t) and Wy (j,t) denote the wavelet coefficients obtained from wavelet transform of X and Y, respectively, at scale j and time t.

The global wavelet coherence at a certain scale j is defined as the time-averaged value of the wavelet coefficients at the scale with the COI. It is estimated by

(21) R 2 j = 1 n t = t 1 t 2 R 2 ( j , t ) ,

where nj is the number of points with COI and nj=t2-t1+1.

Figure 4Spatial distribution/extent of the seven regions, or communities, with similar heavy precipitation event characteristics across India. Black dots indicate the representative grid cells for each of the community identified using the ZP space approach. Terrain characteristics of the Indian subcontinent are shown using the SRTM DEM (in the background).

Global wavelet coherence is a useful measure to examine the common characteristic periodicities between x and y. Grinsted et al. showed the applicability of WC analysis of the association of precipitation with climate variables (Grinsted et al., 2004). A more detailed description of wavelet coherence analysis can be found in Grinsted et al. (2004).

It is important to note that WCA uses the complete, continuous time series for quantifying the linkages between precipitation and climate patterns, whereas MSES first derives extreme events at the different timescales and then uses the synchronization between these events to identify the linkages.

Figure 5Multiscale event synchronization (MSES) between precipitation and climate indices. From top to bottom: Niño 3.4, IOD, NAO, PDO, and AMO. From left to right: community 1 to community 7. MSES values are shown as solid lines, and significant connections (at the 95 % significance level) are marked in grey.


Figure 6Global wavelet coherence (GWC) between precipitation and climate indices. Top to bottom: Niño 3.4, IOD, NAO, PDO, and AMO. From left to right: community 1 to community 7. WCA values are shown as solid lines, and significant connections (at the 95 % significance level) are marked in grey.


4 Results and discussion

4.1 Homogeneous regions and representative grid cells

To reduce the number of pairs of precipitation and climate index time series for finding synchronization, we pool precipitation grid cells with similar heavy precipitation event characteristics into homogeneous regions. These regions and their main physical characteristics are given in Fig. 4. A more detailed discussion of these regions is provided in a previous study (Agarwal et al., 2018a). For each community (C1 to C7), we identify a representative grid cell (black dots in Fig. 4) using the ZP space approach. C1 and C2 (Fig. 4) both are in southern India but are differentiated by topological (elevation, land, coastline and climate regimes) features. C3 has moderate elevation, equatorial grasslands and semi-arid climate regimes. C4 covers almost all of the greenest and most mountainous regions of India (north-eastern India). C5 in north-western India covers dry and lowland areas. C6 in the western coastline is near to both coastlines and low-lying areas with two different climate regimes (arid and humid). C7 is very high mountainous region with alpine climate regimes. Next, we investigate the non-linear linkages between the precipitation time series of the representative cells and the climate indices.

4.2 Linkages between precipitation and climatic patterns at multiple timescales

Figures 5a–e and 6a–e show the MSES values and WCA values between precipitation and the climate indices, respectively. They are given for the five chosen climate indices and extreme precipitation in each of the representative grid cells of the seven homogeneous regions.

Figure 7Schematic map of spatial diversity of Indian precipitation teleconnections at different timescales. (a) Niño 3.4, (b) IOD, (c) NAO, (d) PDO, and (e) AMO. Colours are consistent with the community shown in Fig. 4. Presence of colour (irrespective of magnitude of synchronization) in the community segment indicates significant synchronization between teleconnection and Indian precipitation. Every single segment of circle shows the temporal scale. The cardinal direction has been projected in the background of each circle.


Figure 5a shows a significant association between El Niño–Southern Oscillation (ENSO) and precipitation in all regions of India at the interannual scale. Its strength varies in space and with temporal scale. It is stronger for the south-eastern peninsular (C1, C2, C3 and C4) and decreases notably in the north-western Himalayan (C5, C6 and C7) regions. In the south-eastern peninsula, the highest synchronization for the low (C1), mild (C2) and moderate (C3) elevation regions occurs at the 4-year scale and at the 2-year scale for the high-elevation (C4) region. For the south-eastern regions of India, we observe a significant synchronization at the decadal scale (8–16 years) which is counterintuitive given the interannual timescale of ENSO (D'Arrigo, 2005; McGregor et al., 2013). The analysis based on WCA (Fig. 6a) shows substantially less correlation between precipitation and ENSO in all regions.

Overall, the association between ENSO and precipitation at the interannual scale is coherent with the general understanding that extreme precipitation in India is associated with ENSO (Rajeevan and Pai, 2007). Additionally, our analysis reveals the important spatial variation of this linkage across India, which has not yet been reported before. We find stronger linkages for the regions close to the ocean (south-eastern peninsular comprising C1 to C4) compared to the inland regions with higher elevation (north-western India comprising C5 and C7). The results mentioned above are in congruence with the findings by Guhathakurta et al. (2017) and Mishra et al. (2012). The spatial heterogeneity in the strength of the relationship between ENSO and precipitation may be a result of the tropical convection during the ENSO events (Bansod, 2011). Other studies have confirmed that there is a decrease in the strength of the relationship between precipitation and ENSO events with distance from the ocean. A similar pattern is observed in Mexico where the Niño 3.4 teleconnection is weaker, if not opposite in sign, in northern versus southern Mexico (Hu and Feng, 2002). This observation leads us to the understanding that the ENSO teleconnection is strong in regions of climatologically strong convection.

Interestingly, an association between ENSO and precipitation at the decadal scale has not been reported for India so far. This association might be a consequence of the interdependencies between ENSO and IOD at the decadal scale (Luo et al., 2010). Recently, Izumo et al. (2010), demonstrated that IOD events tend to not only co-occur with ENSO events, but also to lead them through tropospheric biennial oscillation (Pillai and Mohankumar, 2010). MSES has the potential to capture such interdependencies when applied directly to such indices. However, this is beyond the scope of the study.

The synchronization and coherence between the Indian Ocean Dipole (IOD) and precipitation are given in Figs. 5b and 6b, respectively. The non-linear dependence measure points to a significant synchronization at timescales of 8–16 years in the south-eastern regions C1–C4. The rest of the country seems to be unaffected by IOD. The WCA analysis obtains a similar spatial pattern; however, the significant associations occur at shorter timescales (1–4 years, Fig. 6b). Interestingly, with both methods, we cannot find any coupling with the Himalayan region (C7).

The results obtained by MSES and WCA are in accordance with the general understanding that IOD plays a vital role in the Indian monsoon system in the south-eastern regions, i.e. in close proximity to the Indian Ocean, at interannual and decadal scales (Krishnan and Swapna, 2009). This result can be explained by the fact that two of the general conditions for Indian precipitation are the Tropical Easterly Jet and Tropical Westerly Jet (Rai and Dimri, 2017). In the case of occurrence of IOD, the pressure dipole generated between the Tibetan Plateau and Madagascar either strengthens the south-eastern Indian monsoon (positive IOD) or weakens it (negative IOD) (Jiang and Ting, 2017). However, the reason for the association at the decadal scale is not apparent and needs further investigation.

Unlike IOD, NAO demonstrates significant synchronization with precipitation across the entire subcontinent (Fig. 5c). The linkages to the northern regions C4, C5 and C7 are strong and significant at interannual and decadal scales, whereas the southern regions C1, C2, C3 and C6 show weaker linkages. Overall, the strength of synchronization between NAO and Indian extreme precipitation is higher at the decadal scale than at the interannual scales. The comparison of the results obtained by MSES (Fig. 5c) and WCA (Fig. 6c) reveals that the non-linear method shows an increase in the association, particularly in the north-eastern Himalayan foothill region (C4). For some regions, MSES detects linkages which are not found by WCA. For example, in the Himalayan region (C7), MSES shows a significant association at timescales of 4–16 years, whereas WCA shows only a signal just at the significance level at the scale of 16 years. The overall MSES results are in congruence with other studies (Bhatla et al., 2016; Feliks et al., 2013; Goswami et al., 2006), but so far space and scale variation in the associations between NAO and Indian precipitation has not gained attention. The linkages between precipitation and NAO in the northern part of the country might be due to westerly influences from the Eurasian region which are, in turn, strongly affected by NAO. Another explanation (Goswami et al., 2006) suggests that the linkage of NAO and Indian precipitation at higher scales (decadal and beyond) in the northern part of India results from the interdependency of NAO and AMO.

In the case of PDO, we infer a robust decadal synchronization across the entire subcontinent (Fig. 5d). The strength of synchronization varies in space and reaches values of around 0.7 for several regions. By contrast, WCA (Fig. 6d) does not reveal significant associations at the decadal scale except for the eastern coastline (C1) and Himalayan foothills (C4), where values at the boundary with significance are found.

The MSES results agree with Krishnan and Sugi (2003), who demonstrate a strong relationship between PDO and precipitation across the country. The interannual synchronization might be an indirect influence because of the interdependency of PDO and ENSO (Krishnan and Sugi, 2003; Rathinasamy et al., 2014).

The highest strength of synchronization between AMO and Indian precipitation is observed in the north-western and central regions C3 to C6 (Fig. 5e). Weaker associations are detected in the south (C1, C2), whereas no significant synchronization is found for the Himalayan region (C7). The linkages are most prominent at the decadal scale; in some regions significant synchronization at interannual scales is also found. In contrast, WCA shows only weak linkages (Fig. 6e).

Our MSES results confirm the assertion made by Zhang and Delworth (2005), who found an in-phase relationship between Indian precipitation and AMO. A study by Goswami et al. (2006) also unravelled a link between AMO and multidecadal variability of Indian precipitation. However, our study is the first to observe that the strength of the coupling between AMO and precipitation varies according to the different climate regions and is strongest at the decadal scale.

In summary, our findings re-confirm known physics-based associations, thus implicitly affirming the validity of our approach, but also provide new insights into Indian precipitation teleconnections. We find substantial spatial variation in the significant linkages across India and for different timescales (Fig. 7). MSES reveals an appreciable increase in the association between climate patterns and precipitation in most regions when compared to WCA. In some regions, the synchronization values increase by 40 %–50 %. The much higher skill of MSES in detecting associations suggests the presence of non-linear and threshold relationships which cannot be captured by WCA, which is limited to linear processes.

5 Conclusions

A novel non-linear, multiscale approach (MSES) based on wavelets and event synchronization is used for unravelling teleconnection influences on the Indian climate network at multiple timescales. The analysis considers those climate patterns with the highest relevance for Indian precipitation. To understand the spatial heterogeneity, India is sub-divided into homogeneous regions using complex networks. The comparison with wavelet coherence analysis (WCA), the state-of-the-art method in understanding linkages at different timescales, shows a much higher skill for MSES in detecting linkages between climate indices and precipitation. This suggests that there are significant non-linear linkages which are not well captured by linear approaches such as WCA.

The application of MSES to the homogeneous regions, obtained using a complex network approach, allows the spatial diversity in the teleconnection patterns over India to be unravelled. ENSO has a strong influence on precipitation in the south-eastern parts of the country. These regions are also affected by IOD; however, the IOD influence is much weaker compared to ENSO. NAO has a strong connection to extreme precipitation, particularly in the northern regions. The effect of PDO stretches across the whole country, whereas AMO influences precipitation particularly in the arid and semi-arid regions. The substantial variation in precipitation teleconnections across India and across timescales that is unravelled by the proposed method provides an exciting perspective for rainfall forecasting for India and for making better sense of its weather.

Data availability

The authors used a high-resolution (0.25×0.25) monthly gridded precipitation dataset developed by the Indian Meteorological Department (IMD) for the spatial domain of 66.5 to 100 E and 6.5 to 38.5 N, covering the mainland region of India (Pai et al., 2014). The gridded data have been generated from the observations of 6995 gauging stations across India (Pai et al., 2014). Details about these data can be obtained from (last access: 13 August 2019) (homepage/rainfall information).

Author contributions

AA devised the project, the main conceptual ideas, and the proof outline. AA and RM jointly developed the theoretical framework. AA took the lead, implemented the method on the real dataset and performed the analysis and tests. AA arranged and preprocessed the dataset, prepared all the figures and wrote the main text. BM closely supervised the work, helped extensively in rewriting the text for journal submission and helped to improve the figures. Intense discussion with NM and JK helped to reach to appropriate conclusion relevant for the climatological interpretation. RK provided the additional help to draw a meaningful conclusion in terms of climatological interpretations. RS and LC helped me in rewriting the text, climatological interpretations and conceptualizing figures.

Competing interests

The authors declare that they have no conflict of interest.

Special issue statement

This article is part of the special issue “Centennial issue on nonlinear geophysics: accomplishments of the past, challenges of the future”. It is not associated with a conference.


This research was funded by the Deutsche Forschungsgemeinschaft (DFG) (GRK 2043/1) within the graduate research training group “Natural risk in a changing world (NatRiskChange)” at the University of Potsdam (, last access: 13 August 2019). The authors gratefully thank Roopam Shukla and Bedartha Goswami for the helpful suggestion.

Financial support

This research has been supported by the Deutsche Forschungsgemeinschaft (grant no. 2043/1).

Review statement

This paper was edited by Stéphane Vannitsem and reviewed by two anonymous referees.


Abid, M. A., Almazroui, M., Kucharski, F., O'Brien, E., and Yousef, A. E.: ENSO relationship to summer rainfall variability and its potential predictability over Arabian Peninsula region, npj Climate and Atmospheric Science, 1, 20171,, 2018. 

Agarwal, A.: Unraveling spatio-temporal climatic patterns via multi-scale complex networks, Universität Potsdam, 2019. 

Agarwal, A., Maheswaran, R., Kurths, J., and Khosa, R.: Wavelet Spectrum and Self-Organizing Maps-Based Approach for Hydrologic Regionalization – a Case Study in the Western United States, Water Resour. Manag., 30, 4399–4413,, 2016. 

Agarwal, A., Marwan, N., Rathinasamy, M., Merz, B., and Kurths, J.: Multi-scale event synchronization analysis for unravelling climate processes: a wavelet-based approach, Nonlin. Processes Geophys., 24, 599–611,, 2017. 

Agarwal, A., Marwan, N., Maheswaran, R., Merz, B., and Kurths, J.: Quantifying the roles of single stations within homogeneous regions using complex network analysis, J. Hydrol., 563, 802–810,, 2018a. 

Agarwal, A., Maheswaran, R., Marwan, N., Caesar, L., and Kurths, J.: Wavelet-based multiscale similarity measure for complex networks, Eur. Phys. J. B, 91, 296,, 2018b. 

Agarwal, A., Caesar, L., Marwan, N., Maheswaran, R., Merz, B., and Kurths, J.: Network-based identification and characterization of teleconnections on different scales, Sci. Rep., 9, 8808,, 2019a. 

Agarwal, A., Marwan, N., Ozturk, U., and Maheswaran, R.: Unfolding Community Structure in Rainfall Network of Germany Using Complex Network-Based Approach, in: Water Resources and Environmental Engineering II, edited by: Rathinasamy, M., Chandramouli, S., Phanindra, K. B. V. N., and Mahesh, U., 179–193, Springer Singapore, Singapore, 2019b. 

Araghi, A., Mousavi-Baygi, M., Adamowski, J., and Martinez, C.: Association between three prominent climatic teleconnections and precipitation in Iran using wavelet coherence, Int. J. Climatol., 37, 2809–2830,, 2017. 

Ashok, K., Guan, Z., and Yamagata, T.: Impact of the Indian Ocean dipole on the relationship between the Indian monsoon rainfall and ENSO, Geophys. Res. Lett., 28, 4499–4502,, 2001. 

Bansod, S. D.: Interannual variability of convective activity over the tropical Indian Ocean during the El Niño/La Niña events, Int. J. Remote Sens., 32, 5565–5582,, 2011. 

Behera, S. K., Krishnan, R., and Yamagata, T.: Unusual ocean-atmosphere conditions in the tropical Indian Ocean during 1994, Geophys. Res. Lett., 26, 3001–3004,, 1999. 

Bharath, R. and Srinivas, V. V.: Delineation of homogeneous hydrometeorological regions using wavelet-based global fuzzy cluster analysis, Int. J. Climatol., 35, 4707–4727,, 2015. 

Bhatla, R., Singh, A. K., Mandal, B., Ghosh, S., Pandey, S. N., and Sarkar, A.: Influence of North Atlantic Oscillation on Indian Summer Monsoon Rainfall in Relation to Quasi-Binneal Oscillation, Pure Appl. Geophys., 173, 2959–2970,, 2016. 

Blondel, V. D., Guillaume, J.-L., Lambiotte, R., and Lefebvre, E.: Fast unfolding of communities in large networks, J. Stat. Mech.-Theory E., 2008, P10008,, 2008. 

Boers, N., Goswami, B., Rheinwalt, A., Bookhagen, B., Hoskins, B., and Kurths, J.: Complex networks reveal global pattern of extreme-rainfall teleconnections, Nature, 566, 373–377,, 2019. 

D'Arrigo, R.: On the variability of ENSO over the past six centuries, Geophys. Res. Lett., 32,, 2005. 

Ding, R., Kang, S., Vargas, R., Zhang, Y., and Hao, X.: Multiscale spectral analysis of temporal variability in evapotranspiration over irrigated cropland in an arid region, Agr. Water Manag., 130, 79–89,, 2013. 

Dong, X.: Influences of the Pacific Decadal Oscillation on the East Asian Summer Monsoon in non-ENSO years: Influences of the Pacific Decadal Oscillation on the East Asian Summer Monsoon, Atmos. Sci. Lett., 17, 115–120,, 2016. 

Donges, J. F., Zou, Y., Marwan, N., and Kurths, J.: Complex networks in climate dynamics: Comparing linear and nonlinear network construction methods, The European Physical Journal Special Topics, 174, 157–179,, 2009. 

Ekhtiari, N., Agarwal, A., Marwan, N., and Donner, R. V.: Disentangling the multi-scale effects of sea-surface temperatures on global precipitation: A coupled networks approach, Chaos: An Interdisciplinary Journal of Nonlinear Science, 29, 063116,, 2019. 

Feliks, Y., Groth, A., Robertson, A. W., and Ghil, M.: Oscillatory Climate Modes in the Indian Monsoon, North Atlantic, and Tropical Pacific, J. Climate, 26, 9528–9544,, 2013. 

Feng, Q. Y., Vasile, R., Segond, M., Gozolchiani, A., Wang, Y., Abel, M., Havlin, S., Bunde, A., and Dijkstra, H. A.: ClimateLearn: A machine-learning approach for climate prediction using network measures, Geosci. Model Dev. Discuss.,, 2016. 

Fortunato, S.: Community detection in graphs, Phys. Rep., 486, 75–174,, 2010. 

Goswami, B. N. and Krishnan, R.: Opportunities and challenges in monsoon prediction in a changing climate, Clim. Dynam., 41, 1–1,, 2013. 

Goswami, B. N., Madhusoodanan, M. S., Neema, C. P., and Sengupta, D.: A physical mechanism for North Atlantic SST influence on the Indian summer monsoon, Geophys. Res. Lett., 33,, 2006. 

Grinsted, A., Moore, J. C., and Jevrejeva, S.: Application of the cross wavelet transform and wavelet coherence to geophysical time series, Nonlin. Processes Geophys., 11, 561–566,, 2004. 

Guhathakurta, P., Menon, P., Inkane, P. M., Krishnan, U., and Sable, S. T.: Trends and variability of meteorological drought over the districts of India using standardized precipitation index, J. Earth Syst. Sci., 126, 120,, 2017. 

Guimerà, R. and Amaral, L. A. N.: Cartography of complex networks: modules and universal roles, J. Stat. Mech.-Theory E., 2005, P02001,, 2005. 

Halverson, M. J. and Fleming, S. W.: Complex network theory, streamflow, and hydrometric monitoring system design, Hydrol. Earth Syst. Sci., 19, 3301–3318,, 2015. 

Hannachi, A., Jolliffe, I. T., and Stephenson, D. B.: Empirical orthogonal functions and related techniques in atmospheric science: A review, Int. J. Climatol., 27, 1119–1152,, 2007. 

Harenberg, S., Bello, G., Gjeltema, L., Ranshous, S., Harlalka, J., Seay, R., Padmanabhan, K., and Samatova, N.: Community detection in large-scale networks: a survey and empirical evaluation: Community detection in large-scale networks, Wiley Interdisciplinary Reviews: Computational Statistics, 6, 426–439,, 2014. 

Hu, Q. and Feng, S.: Interannual Rainfall Variations in the North American Summer Monsoon Region: 1900–98*, J. Climate, 15, 1189–1202,<1189:IRVITN>2.0.CO;2, 2002. 

Hu, W. and Si, B. C.: Technical note: Multiple wavelet coherence for untangling scale-specific and localized multivariate relationships in geosciences, Hydrol. Earth Syst. Sci., 20, 3183–3191,, 2016. 

Izumo, T., Vialard, J., Lengaigne, M., de Boyer Montegut, C., Behera, S.K., Luo, J.-J., Cravatte, S., Masson, S., and Yamagata, T.: Influence of the state of the Indian Ocean Dipole on the following year’s El Niño, Nat. Geosci., 3, 168–172,, 2010. 

Jiang, X. and Ting, M.: A Dipole Pattern of Summertime Rainfall across the Indian Subcontinent and the Tibetan Plateau, J. Climate, 30, 9607–9620,, 2017. 

Katul, G., Lai, C.-T., Schäfer, K., Vidakovic, B., Albertson, J., Ellsworth, D., and Oren, R.: Multiscale analysis of vegetation surface fluxes: from seconds to years, Adv. Water Resour., 24, 1119–1132,, 2011. 

Krishnamurthy, L. and Krishnamurthy, V.: Teleconnections of Indian monsoon rainfall with AMO and Atlantic tripole, Clim. Dynam., 46, 2269–2285,, 2016. 

Krishnan, R. and Sugi, M.: Pacific decadal oscillation and variability of the Indian summer monsoon rainfall, Clim. Dynam., 21, 233–242,, 2003. 

Krishnan, R. and Swapna, P.: Significant Influence of the Boreal Summer Monsoon Flow on the Indian Ocean Response during Dipole Events, J. Climate, 22, 5611–5634,, 2009. 

Kumar, K. K., Rajagopalan, B., Hoerling, M., Bates, G., and Cane, M.: Unraveling the Mystery of Indian Monsoon Failure During El Nino, Science, 314, 115–119,, 2006. 

Lancichinetti, A. and Fortunato, S.: Community detection algorithms: A comparative analysis, Phys. Rev. E, 80, 056117,, 2009. 

Luo, J.-J., Zhang, R., Behera, S. K., Masumoto, Y., Jin, F.-F., Lukas, R., and Yamagata, T.: Interaction between El Niño and Extreme Indian Ocean Dipole, J. Climate, 23, 726–742,, 2010. 

Luterbacher, J., Xoplaki, E., Casty, C., Wanner, H., Pauling, A., Küttel, M., Rutishauser, T., Brönnimann, S., Fischer, E., Fleitmann, D., Gonzalez-Rouco, F. J., García-Herrera, R., Barriendos, M., Rodrigo, F., Gonzalez-Hidalgo, J. C., Saz, M. A., Gimeno, L., Ribera, P., Brunet, M., Paeth, H., Rimbu, N., Felis, T., Jacobeit, J., Dünkeloh, A., Zorita, E., Guiot, J., Türkes, M., Alcoforado, M. J., Trigo, R., Wheeler, D., Tett, S., Mann, M. E., Touchan, R., Shindell, D. T., Silenzi, S., Montagna, P., Camuffo, D., Mariotti, A., Nanni, T., Brunetti, M., Maugeri, M., Zerefos, C., Zolt, S. D., Lionello, P., Nunes, M. F., Rath, V., Beltrami, H., Garnier, E., and Ladurie, E. L. R.: Chapter 1 Mediterranean climate variability over the last centuries: A review, in: Developments in Earth and Environmental Sciences, 4, 27–148, Elsevier, 2006. 

McGregor, S., Timmermann, A., England, M. H., Elison Timm, O., and Wittenberg, A. T.: Inferred changes in El Niño–Southern Oscillation variance over the past six centuries, Clim. Past, 9, 2269–2284,, 2013. 

Miralles, D. G., Teuling, A. J., van Heerwaarden, C. C., and Vilà-Guerau de Arellano, J.: Mega-heatwave temperatures due to combined soil desiccation and atmospheric heat accumulation, Nat. Geosci., 7, 345–349,, 2014. 

Mishra, V., Smoliak, B. V., Lettenmaier, D. P., and Wallace, J. M.: A prominent pattern of year-to-year variability in Indian Summer Monsoon Rainfall, P. Natl. Acad. Sci. USA, 109, 7213–7217,, 2012. 

Mitra, C., Kurths, J., and Donner, R. V.: Rewiring hierarchical scale-free networks: Influence on synchronizability and topology, EPL-Europhys. Lett., 119, 30002,, 2017. 

Mokhov, I. I., Smirnov, D. A., Nakonechny, P. I., Kozlenko, S. S., and Kurths, J.: Relationship between El-Niño/Southern Oscillation and the Indian monsoon, Izvestiya, Atmos. Ocean. Phys., 48, 47–56,, 2012. 

Newman, M. E. J.: Detecting community structure in networks, Eur. Phys. J. B, 38, 321–330,, 2004. 

Newman, M. E. J.: Modularity and community structure in networks, P. Natl. Acad. Sci. USA, 103, 8577–8582,, 2006. 

Okin, G. S., Parsons, A. J., Wainwright, J., Herrick, J. E., Bestelmeyer, B. T., Peters, D. C., and Fredrickson, E. L.: Do Changes in Connectivity Explain Desertification?, BioScience, 59, 237–244,, 2009. 

Ouachani, R., Bargaoui, Z., and Ouarda, T.: Power of teleconnection patterns on precipitation and streamflow variability of upper Medjerda Basin, Int. J. Climatol., 33, 58–76,, 2013. 

Ozturk, U., Marwan, N., Korup, O., Saito, H., Agarwal, A., Grossman, M. J., Zaiki, M., and Kurths, J.: Complex networks for tracking extreme rainfall during typhoons, Chaos: An Interdisciplinary Journal of Nonlinear Science, 28, 075301,, 2018. 

Pai, D. S., Sridhar, L., Rajeevan, M., Sreejith, O. P., Satbhai, N. S., and Mukhopadyay, B.: Development of a new high spatial resolution (0.25×0.25) Long Period (1901–2010) daily gridded rainfall data set over India and its comparison with existing data sets over the region, Mausam, 65, 1–18, 2014. 

Pai, D. S., Sridhar, L., Badwaik, M. R., and Rajeevan, M.: Analysis of the daily rainfall events over India using a new long period (1901–2010) high resolution (0.25×0.25) gridded rainfall data set, Clim. Dynam., 45, 755–776,, 2015. 

Paluš, M.: Cross-Scale Interactions and Information Transfer, Entropy, 16, 5263–5289,, 2014. 

Percival, D. B.: Analysis of Geophysical Time Series Using Discrete Wavelet Transforms: An Overview, in: Nonlinear Time Series Analysis in the Geosciences, vol. 112, edited by: Donner, R. V. and Barbosa, S. M., 61–79, Springer Berlin Heidelberg, Berlin, Heidelberg, 2008. 

Percival, D. B and Walden, A. T.: Wavelet methods for time series analysis, 4, Cambridge university press, 2000. 

Peters, D. P. C., Pielke, R. A., Bestelmeyer, B. T., Allen, C. D., Munson-McGee, S., and Havstad, K. M.: Cross-scale interactions, nonlinearities, and forecasting catastrophic events, P. Natl. Acad. Sci. USA, 101, 15130–15135,, 2004. 

Peters, D. P. C., Bestelmeyer, B. T., and Turner, M. G.: Cross–Scale Interactions and Changing Pattern–Process Relationships: Consequences for System Dynamics, Ecosystems, 10, 790–796,, 2007. 

Pillai, P. A. and Mohankumar, K.: Individual and combined influence of El Niño-Southern Oscillation and Indian Ocean Dipole on the Tropospheric Biennial Oscillation, Q. J. Roy. Meteor. Soc., 136, 297–304,, 2010. 

Quiroga, R. Q., Kreuz, T., and Grassberger, P.: Event synchronization: A simple and fast method to measure synchronicity and time delay patterns, Phys. Rev. E, 66, 041904-1,, 2002. 

Rai, P. and Dimri, A. P.: Effect of changing tropical easterly jet, low level jet and quasi-biennial oscillation phases on Indian summer monsoon: TEJ, LLJ and QBO phases and Indian summer monsoon, Atmos. Sci. Lett., 18, 52–59,, 2017. 

Rajeevan, M. and Pai, D. S.: On the El Niño-Indian monsoon predictive relationships, Geophys. Res. Lett., 34, L04704,, 2007. 

Rathinasamy, M., Khosa, R., Adamowski, J., ch, S., Partheepan, G., Anand, J., and Narsimlu, B.: Wavelet-based multiscale performance analysis: An approach to assess and improve hydrological models, Water Resour. Res., 50, 9721–9737,, 2014. 

Rubinov, M. and Sporns, O.: Weight-conserving characterization of complex functional brain networks, NeuroImage, 56, 2068–2079,, 2011.  

Shukla, R., Agarwal, A., Sachdeva, K., Kurths, J., and Joshi, P. K.: Climate change perception: an analysis of climate change and risk perceptions among farmer types of Indian Western Himalayas, Climatic Change, 152, 103–119,, 2018. 

Sivakumar, B., Singh, V. P., Berndtsson, R., and Khan, S. K.: Catchment Classification Framework in Hydrology: Challenges and Directions, J. Hydrol. Eng., 20, A4014002,, 2015. 

Steinhaeuser, K., Ganguly, A. R., and Chawla, N. V.: Multivariate and multiscale dependence in the global climate system revealed through complex networks, Clim. Dynam., 39, 889–895,, 2012. 

Tan, X., Gan, T. Y., and Shao, D.: Wavelet analysis of precipitation extremes over Canadian ecoregions and teleconnections to large-scale climate anomalies: Large Precipitation and Climate Anomalies, J. Geophys. Res.-Atmos., 121, 14469–14486,, 2016. 

Torrence, C. and Compo, G. P.: A Practical Guide to Wavelet Analysis, B. Am. Meteorol. Soc., 79, 61–78,<0061:APGTWA>2.0.CO;2, 1998. 

Tsonis, A. A., Swanson, K. L., and Roebber, P. J.: What Do Networks Have to Do with Climate?, B. Am. Meteorol. Soc., 87, 585–595,, 2006. 

Tsonis, A. A., Wang, G., Swanson, K. L., Rodrigues, F. A., and da Fontura Costa, L.: Community structure and dynamics in climate networks, Clim. Dynam., 37, 933–940,, 2011. 

Xoplaki, E., González-Rouco, J. F., Luterbacher, J., and Wanner, H.: Wet season Mediterranean precipitation variability: influence of large-scale dynamics and trends, Clim. Dynam., 23, 63–78,, 2004. 

Zhang, R. and Delworth, T. L.: Simulated Tropical Response to a Substantial Weakening of the Atlantic Thermohaline Circulation, J. Climate, 18, 1853–1860,, 2005. 

Short summary
We examined the spatial diversity of Indian rainfall teleconnection at different timescales, first by identifying homogeneous communities and later by computing non-linear linkages between the identified communities (spatial regions) and dominant climatic patterns, represented by climatic indices such as El Nino–Southern Oscillation, Indian Ocean Dipole, North Atlantic Oscillation, Pacific Decadal Oscillation and Atlantic Multi-Decadal Oscillation.