Articles | Volume 25, issue 11
Research article
04 Nov 2021
Research article |  | 04 Nov 2021

Spatiotemporal and cross-scale interactions in hydroclimate variability: a case-study in France

Manuel Fossa, Bastien Dieppois, Nicolas Massei, Matthieu Fournier, Benoit Laignel, and Jean-Philippe Vidal

Understanding how water resources vary in response to climate at different temporal and spatial scales is crucial to inform long-term management. Climate change impacts and induced trends may indeed be substantially modulated by low-frequency (multi-year) variations, whose strength varies in time and space, with large consequences for risk forecasting systems. In this study, we present a spatial classification of precipitation, temperature, and discharge variability in France, based on a fuzzy clustering and wavelet spectra of 152 near-natural watersheds between 1958 and 2008. We also explore phase–phase and phase–amplitude causal interactions between timescales of each homogeneous region. A total of three significant timescales of variability are found in precipitation, temperature, and discharge, i.e., 1, 2–4, and 5–8 years. The magnitude of these timescales of variability is, however, not constant over the different regions. For instance, southern regions are markedly different from other regions, with much lower (5–8 years) variability and much larger (2–4 years) variability. Several temporal changes in precipitation, temperature, and discharge variability are identified during the 1980s and 1990s. Notably, in the southern regions of France, we note a decrease in annual temperature variability in the mid 1990s. Investigating cross-scale interactions, our study reveals causal and bi-directional relationships between higher- and lower-frequency variability, which may feature interactions within the coupled land–ocean–atmosphere systems. Interestingly, however, even though time frequency patterns (occurrence and timing of timescales of variability) were similar between regions, cross-scale interactions are far much complex, differ between regions, and are not systematically transferred from climate (precipitation and temperature) to hydrological variability (discharge). Phase–amplitude interactions are indeed absent in discharge variability, although significant phase–amplitude interactions are found in precipitation and temperature. This suggests that watershed characteristics cancel the negative feedback systems found in precipitation and temperature. This study allows for a multi-timescale representation of hydroclimate variability in France and provides unique insight into the complex nonlinear dynamics of this variability and its predictability.

1 Introduction

Hydroclimate variability represents the spatiotemporal evolution of hydrological (e.g., discharge and groundwater level) and climate variables (e.g., precipitation and temperature) which are directly impacting hydrological variability. Studying how hydrological variables react to climate variability and change is a major challenge for society, in particular for water resource management and flood and drought mitigation planning (IPCC2007, 2014, 2021). However, hydrological variability is expressed at multiple timescales (Labat2006; Schaefli et al.2007; Massei et al.2007, 2017), for which the driving mechanisms remain poorly characterized and understood. As suggested in Blöschl et al. (2019), understanding the spatiotemporal scaling, i.e., how the general dynamics driving hydrological variability change at spatial and temporal scales, represents a major challenge toward improved prediction systems (Gentine et al.2012). Understanding spatiotemporal scaling required to identify regions, i.e., the maximum spatial scale in which the dynamics remain unchanged despite its nonlinearity, is critical (Hubert2001). Hydrological variability is, by definition, nonlinear (Labat2000; Lavers et al.2010; McGregor2017) as it results from complex interactions between atmospheric dynamics and catchment properties that may vary at different timescales (e.g., soil characteristics, water table, karstic systems, and vegetation covers; Gudmundsson et al.2011; Sidibe et al.2019). Such interactions between processes at different timescales, i.e., cross-scale interactions (Palus2014; Jajcay et al.2018), have never been studied to further understand hydrological variability. It has also been shown that hydroclimate variability is inherently nonstationary, with the time dependence of the mean and variance due to changes in the controlling factors (e.g., Coulibaly and Burn2004; Labat2006; Dieppois et al.2013, 2016; Massei et al.2017). This results in difficulties in characterizing and predicting the hydrological variability at different spatiotemporal scales (Gentine et al.2012; Blöschl et al.2019).

While different timescales have been identified in hydrological variability (Coulibaly and Burn2004; Labat2006; Dieppois et al.2013, 2016; Massei et al.2017), very little has been done to (i) explore how spatially coherent those timescales are and (ii) identify regions in which the statistical characteristics of all ranges of variability remain unchanged. Studying 231 stream gauges throughout the world, Labat (2006) highlighted different timescales of discharge variability over the different continents. At the regional scale, Smith et al. (1998) established a clustering of 91 USA stream gauges based on their global wavelet spectra, i.e., dominant timescales, and found five homogeneous regions. Similarly, Anctil and Coulibaly (2004) and Coulibaly and Burn (2004) established a clustering of southern Quebec and Canadian streamflow, based on the timing of both the 2–3- and 3–6-year timescales. In Europe, Gudmundsson et al. (2011) identified different regions according to the magnitude of decadal discharge variability. In France, such a clustering, based on time frequency patterns of discharge variability, and its relation to climate variability (e.g., precipitation and temperature), has not yet been explored. In addition, all studies mentioned above either isolated particular timescales of variability or averaged the variability across timescales (e.g., global wavelet spectra), which is equivalent to a linearization of the system (Hubert et al.1989). These studies thus ignore potential feedback mechanisms, e.g., between soil moisture, precipitation, and temperature (Materia et al.2021; Ardilouze et al.2020; Bellucci et al.2015). In the presence of feedback mechanisms, interactions occur over different timescales, which are also called cross-scale interactions (Christophe Bouton2017). However, while studying cross-scale interactions has gained increasing interest in other fields, such as neurosciences (e.g., Onslow et al.2014; Wang et al.2014), cross-scale interactions are poorly understood in climate and hydrological sciences. New strategies have recently been developed to facilitate such studies (Jajcay et al.2018). Cross-scale interactions are, however, very relevant to hydroclimate studies, in particular when searching for climate drivers (or predictors) of hydrological signals, as they will reveal climate timescales that are causality linked to each timescale of hydrological variability.

In this study, we investigate the spatial homogeneity of hydroclimate variability in France across timescales. We aim at identifying homogeneous regions according to specific time frequency patterns. From the determination of homogeneous regions of hydroclimate variability, we will explore cross-scale interactions that may result from feedback processes between catchment properties and hydroclimate variability.

This study, therefore, has major implications for the comprehension of hydroclimate dynamics and their interactions with large-scale climate drivers and catchment properties. In addition, as recently suggested in Scaife and Smith (2018), improved characterization of the different timescales of variability and their interactions could help optimize ensemble-based hydrological forecasting systems through identifying climate ensemble members that better match the observed realization.

The work is divided into the following sections. Data and methods are introduced in Sect. 2. In Sect. 3, we establish homogeneous regions for precipitation, temperature, and discharge variability based on their time frequency patterns and then explore cross-scale interactions for each region of homogeneous variability in precipitation, temperature, and discharge. Finally, discussions of the main results and conclusions are provided in Sect. 4.


Figure 1Location of stream gauges (gray dots), corresponding watersheds (pale red; Brigode et al.2020), hydrographic network (blue lines; Pella et al.2012), and orography in Safran data set (grayscale; Vidal et al.2010)

Figure 2The workflow of this study is repeated for precipitation, temperature, and discharge data sets. (a) A total of 152 near-natural watersheds are selected. (b) Each watershed is represented by a monthly time series from 1968 to 2008, discharge is measured at a gauging station, and precipitation and temperature are averaged over the watershed. (c) The continuous wavelet transform of each time series is computed, representing the timescale-dependent, nonstationary variability in each watershed. (d) The similarity between all 152 continuous wavelet transforms is computed and represented as a distance matrix. (e) Similar watersheds are grouped together into regions of homogeneous variability using a fuzzy clustering algorithm (top). For each region, the relative importance of each timescale-dependent variability is represented with a global wavelet spectrum, and global wavelet spectra are superimposed for all regions (middle). The continuous wavelet spectra of the regions are superimposed (bottom). For clarity, only timescale and time locations with the most significant variability are shown (colored circles). (f) For each region, cross-timescale interactions are computed. The phase–phase interactions (top) identify any timescale's phase that conditions another timescale's phase. The phase–amplitude interactions (bottom) characterize any timescale's phase that conditions the amplitude of another timescale.

2 Data and methodology

2.1 Hydrological and climate data

The data consist of precipitation, temperature, and discharge time series located over 152 watersheds (Figs. 1, 2a and b). Discharge time series were extracted from French reference hydrometric network compiled by Giuntoli et al. (2013). This network of stations identifies near-natural watersheds (i.e., with negligible anthropogenic modifications) with long-term, high-quality hydrometric data. According to Giuntoli et al. (2013), this subset of stations does not show abrupt changes and trends that could have resulted from anthropogenic influence. The period of 1968–2008 was chosen by Giuntoli et al. (2013) as being the best trade-off in terms of data availability over the different regions. Here, this database was further subset to 152 watersheds in order to select complete monthly time series (i.e., without missing values) only (Fig. 1). Precipitation and temperature data have been estimated from the 8 km grid Safran surface reanalysis data set (Vidal et al.2010) and have been subset to a common period (1968–2008). For this study, precipitation and temperature have been averaged over each watershed area (Caillouet et al.2017). Each station is thus representative of one watershed.

2.2 Methods

The methodology described below (and summarized in the workflow in Fig. 2) is applied to precipitation, temperature, and discharge data sets.

2.2.1 Continuous wavelet transforms

For each of the 152 watersheds, continuous wavelet analysis is used to identify at which timescales and time locations the amplitude of variability (i.e., local variance) is the strongest (Fig. 2c; Torrence and Compo1998). Here we employ the words timescale and frequency interchangeably, though frequency implies a periodic variability, which is not a necessary condition in continuous wavelet analysis. For any finite energy signal x, it is possible to obtain a time frequency representation by projecting the time series onto a function called the mother wavelet, which quantifies the amplitude of the time series variability at a given timescale and time location. This mother wavelet can be translated in time to quantify the variability at precise time locations but can also be scaled so that variability at different timescales can be quantified as well (Torrence and Compo1998; Grinsted et al.2004). A mother wavelet at a timescale a and time location b is called a daughter wavelet. Daughter wavelets are calculated as follows:

(1) ψ a , b = 1 a ψ t - b a .

The left-hand side (LHS) term is the daughter wavelet of scale a and time translation b at time t. For the sake of simplicity, we will refer to b as the time location. The first right-hand side (RHS) term is the scaling of the mother wavelet ψ, and the last one is the time translation. The projection of the signal onto each scale a takes the following form:

(2) WT ψ [ x ] ( a , b ) = x , ψ a , b = R x ( t ) ψ a , b ( t ) d t .

The LHS term contains the wavelet coefficients, WT, i.e., how large the amplitude of variability at the timescale a and time location b is. If the mother wavelet (and, hence, the daughter wavelets as well) is complex, wavelet coefficients are complex as well, and both the amplitude and instantaneous phase of the time series can be computed around time location b and timescale a. Wavelet coefficients represent the inner product of the signal, the daughter wavelet of scale a, and the time location b (center). The norm of their square is called the wavelet power and represents the amplitude of the oscillation of signal x at scale a and centered on time location b. As it is impossible to capture the best resolution in both timescale and time location simultaneously, here we used a Morlet mother wavelet (of the order of 6), which offers a good trade-off between the detection of scales and localization of the oscillations in time (Torrence and Compo1998). Visualization of a continuous wavelet transform is called a scalogram. Figure 2c shows a collection of scalograms, with the time location on the horizontal axis and timescale on the vertical axis. Yellow colors show the timescales and time locations when the amplitude of the time series' variability is maximum. A major advantage of continuous wavelet transform, compared to other signal analysis methods such as the Fourier transform, is that wavelet analysis takes nonstationarity into account. Nonstationarity is the time location dependence of both the mean and variance of a time series.

Because the daughter wavelet translates and scales up, overlaps in time and frequency can occur, and wavelet coefficients can be overestimated, requiring statistical significance tests (Torrence and Compo1998). This redundancy may give rise to peaks in the wavelet coefficients (meaning that high variability is detected) even in the case of a random noise (Ge2007). Torrence and Compo (1998) used Monte Carlo simulations to assess the statistical significance of the continuous wavelet transform of their time series. Just as with any other statistical analysis, the performance of statistical tests is an open debate. It has been shown that, while the significance test of Torrence and Compo (1998) may not unravel all significant wavelet coefficients, their false detection rate is low for as long as the mother wavelet chosen is adapted to the time series (Ge2007). By using a Morlet wavelet of the order of 6, we ensure that such statistical significance tests keep the false detection rate low (Ge2007). Because the significance test's aim is to ensure that large peaks are exceeding the range of variability that would occur in a random noise, statistically significant wavelet coefficients are always those with large wavelet coefficients values for any given timescale.

For the reminder of this study, the terms intraseasonal, annual, and interannual refer to variations at <1-, 1-, 2–4-, and 5–8-year timescales, respectively.

2.2.2 Image Euclidean distance clustering

After each watershed wavelet spectrum is computed, we estimate the similarity between them, i.e., how similar the variability is, for given scales and time locations, among all wavelet spectra (Fig. 2d). Similarities between wavelet spectra are estimated from the entire wavelet spectrum, and not only on statistically significant signals, to guarantee more consistent comparison between spectra. Distances between two-dimensional data, such as wavelet spectra, are estimated using Euclidean distance between pairwise points (pED; i.e., computing f2(xi,yi)-f1(xi,yi)). However, such a procedure has no neighborhood notion, making it impossible to account for globally similar shapes (Wang et al.2005). To avoid this issue, we used the image Euclidean distance calculation method (hereinafter IEDC) developed by Wang et al. (2005). The IEDC method modifies the pED equation in the following two ways (Wang et al.2005): (i) the distance between pixel values is computed not only pairwise but for all indices, and (ii) a Gaussian filter, a function of the spatial distance between pixels, is applied. The Gaussian filter then applies less weight to the computed distance between very close and far apart pixels, while emphasizing on medium-spaced ones (Wang et al.2005).

2.2.3 Fuzzy clustering

Fuzzy clustering has then been used to cluster the different watersheds based on their similarities (Fig. 2e). Fuzzy clustering is a soft clustering method (Dunn1973). While soft clustering spreads membership over all clusters with varying probability, hard clustering attributes each station one, and only one, cluster membership. Soft clustering is therefore better suited to situations when the spatial variability, originating from different stations' hydroclimate characteristics, is smooth. For instance, precipitation and temperature patterns are unlikely to change suddenly from one station to a neighboring one and in turn, be markedly different from the next neighbor (Moron et al.2007; Hannaford et al.2009; Rahiz and New2012). As such, several stations tend to show transitional or hybrid patterns and can potentially be members of different clusters, limiting the robustness of hard clustering procedure (Liu and Graham2018).

Fuzzy clustering performance is determined by the ability of the algorithms to recognize hybrid stations (i.e., stations incorporating multiple features from different patterns observed in other coherent regions), while allowing for a clear determination of the membership of stations with unique features (Kaufman and Rousseeuw1990). Here, we used the FANNY algorithm (Kaufman and Rousseeuw1990), which has been shown to be flexible and to offer the possibility to adapt the clustering to the data with optimal performance (Liu and Graham2018). In addition, rather than setting the number of clusters arbitrarily, we used an estimation of the optimum number of clusters by first computing a hard clustering method, namely consensus clustering (Monti et al.2003). Thus, the number of clusters providing the best stability (i.e., the minimal changes in membership when adding new individuals) is considered optimal, as recommended in Şenbabaoǧlu et al. (2014). The different clusters' memberships are then mapped to discuss the spatial coherence of each hydroclimate variable.

2.2.4 Cross-scale interactions

For each variable and each cluster, cross-scale interactions are explored (Fig. 2f). Cross-scale interactions refer to phase–phase and phase–amplitude couplings between timescales of a given time series (Paluš2014; Scheffer-Teixeira and Tort2016). Here, coupling means that the state (either phase or amplitude) of a signal y is dependent on the state of a signal x and describes causal relationship (Granger1969; Pikovsky et al.2001), which refers to the information transfer from a timescale of one signal to another.

Figure 3A system with directional cross-scale interactions. (a) A variable f(t) made of two components, X and Y, connected through CXY and CYX in a perturbation–dampening scheme, so that f(t)=X(t)-Y(t). Both X and Y receive inputs ϕX and ϕY, respectively. CXX allows X to grow first. Depending on both inputs and connections, some phase–phase or phase–amplitude interactions between X and Y can occur. (b) An example of a phase–phase interaction, with every fourth ridge of YPP coinciding with a ridge of X, with fPP(t)=X(t)-YPP(t) (top, middle, and bottom panels, respectively). (c) An example of a phase–amplitude interaction. X and YPA only interact when X reaches a ridge, in which case the YPA amplitude, if lowered, yields fPA(t) (top, middle, and bottom panels, respectively; adapted from Onslow et al.2014).


Figure 3 describes the necessary setting and characteristics of cross-scale interactions. A variable (f) measures the dynamics of a system (e.g., precipitation or temperature variability). This system is modeled as a coupling of two components (X and Y). The components interact with each other in a perturbation–dampening (X and Y, respectively) scheme so that f(t)=X(..,t)-Y(..,t) (Fig. 3a). The interactions between the components occur through the connections CXY and CYX, with a given strength (here C..=2), and this perturbation–dampening interaction forms a negative feedback (i.e., an increase in X activity triggers Y activity and dampens X activity, which, in turn, lowers X activity, thus lowering Y activity, thus allowing X activity to increase again, and so on; Fig. 3a). The connection CXX enables X to grow first before Y dampens it. This connection CXX forms a positive feedback, i.e., increases in X activity will be more severe as X activity is high. Both X and Y receive inputs ϕX and ϕY from driving processes (e.g., moisture advection and convective processes; Fig. 3a). Depending on both the mean and timescales of ϕX and ϕY, the strength of CXX, CXY and CYX, and X and Y may show coupled behaviors. For instance, in Fig. 3b, every fourth ridge of YPP(t) is synchronized with a ridge of X(t) (Fig. 3b; top and middle panels), thus forming a phase–phase interaction. The direction of the interactions depends on the inputs (ϕX, ϕY) and the connections (CXX, CXY, and CYX). fPP(t) is the difference between X(t) and YPP (Fig. 3b; bottom panel). Because the interaction between X and Y depends on both inputs and connections, interactions may lead to a cross-scale relationship only for certain values of X or Y (Fig. 3c). Thus, depending on the phase of either X or Y, the amplitude of the driven component may increase/decrease when the cross-scale interaction takes place and return to normal when it is out of phase compared to the driving component. This describes a phase–amplitude interaction. In Fig. 3b, YPA amplitude decreases when X is at its maximum (i.e., when its phase is a ridge; Fig. 3c; top and middle panels). Similar to the phase–phase interaction, fPA(t) is the difference between X(t) and YPA(t) (Fig. 3c; bottom panel). Because phase and amplitude are very dependent on inputs ϕX and ϕY, connections between spatially distant physical processes are likely to give rise to phase–amplitude interactions (Nandi et al.2019). In summary, there are three elements needed for cross-scale interactions between multiple, coupled processes, to arise, namely that (i) an oscillating forcing ϕX must drive X and an additional forcing ϕY on Y may also be present, (ii) X must have positive feedback on itself so that it grows faster than Y;, and(iii) Y must show a dampening effect on X (XY or YX negative feedback). The presence and characteristics of cross-scale interactions depend on the strength and frequencies of ϕX,ϕY, intrinsic frequencies of X,Y, and coupling strengths CXY, CYX, CXX, and CYY (Fig. 3). Thus, the detection of cross-scale interactions in time series is an indication of the presence of all those characteristics in the hydroclimate system, e.g., precipitation–land processes, which helps in investigating potential processes at play.

The balance between X and Y determines if the feedback is either positive or negative (Peters et al.2007).

Note that cross-scale interactions can occur from large-scale to small-scale processes and vice-versa. For instance, atmospheric circulation at seasonal timescales influences interannual and decadal timescales, which, in turn, influences seasonal variations (Hannachi et al.2017).

Following Palus (2014) and Jajcay et al. (2018), we chose the conditional mutual information (CMI) surrogates method, combined with wavelet transforms. First, using a Morlet mother wavelet, the instantaneous phase and amplitude at time t and scale s of the signal are obtained. Next, the conditional mutual information, I(ϕx(t);ϕy(t+τ)-ϕy(t)|ϕy(t)), for the phase, and I(ϕx(t);Ay(t+τ)|Ay(t)),Ay(t-η)),Ay(t-2η)), for the amplitude, is computed. In the case of phase–phase relationships, the CMI measures how much the present phase of x contains information about the future phase of y, knowing the present value of y. Phase–phase interactions can be uni- or bi-directional. It is possible for a single timescale to drive another, which, in turn, drives back the original one, describing feedback interactions. For phase–amplitude relationships, CMI measures how much the present phase of x contains information of the future amplitude of y, knowing the present and past values of y. The statistical significance of the CMI measure is assessed using 5000 phase-randomized surrogates having the same Fourier spectrum and mean and standard deviation as the original time series, as seen in Ebisuzaki (1997). Palus (2014) has shown that this number of surrogates is ideal for statistical significance in the context of hydroclimate time series. The computational cost is, however, high, with approximately 1 week of computing for a time series of 50 years on a 32 core Xeon computer. The present computations were done on the Myria cluster, hosted by the Centre Régional Informatique et d'Applications Numérique de Normandie (; last access: 7 January 2021).

3 Spatiotemporal clustering of hydrological variability

The wavelet transform corresponding to each watershed's monthly time series have been computed, and all 152 watersheds' wavelet transforms have been checked for similarities using IEDC fuzzy clustering to identify and characterize homogeneous regions of hydroclimate variability over France. Once the homogeneous regions had been identified, an average time series for each region was computed. The global wavelet spectrum of this time series quantified the total variance expressed at each timescale, while its wavelet spectrum characterized how this variance is distributed in the time (location) and frequency (scale) domains. In addition, so as to focus on interannual timescales, we computed the wavelet spectrum of the time series filtered at the annual time step. Cross-scale interactions were then investigated for each homogeneous region.

Figure 4Clustering of precipitation time frequency variability in France. (a) Classification map of the watersheds. Pie chart slices show the three highest probability memberships. Pie charts denote fuzzy clustering memberships. (b) Global wavelet spectra of homogeneous regions. (c) Wavelet spectra of homogeneous regions. For clarity, only timescales and time locations that are 95 % statistically significant and with the largest variability are shown (colored circles).

3.1 Precipitation

3.1.1 Time frequency patterns

The following seven regions with homogeneous time frequency patterns are identified (Fig. 4a): northwestern (green), northeastern (blue), central north (red), central western (pink), central eastern (black), southwestern (yellow), and southeastern (dark green). Figure 4a shows that all watersheds converge toward singular clusters, meaning that all regions are highly coherent (i.e., pie charts in Fig. 4a show one dominant color).

In all regions, precipitation is varying at different timescales, ranging from intraseasonal to interannual scales (i.e., 2–8 years; Fig. 4b). Southwestern and southeastern regions are dominated by annual (1 year) variability, while their interannual variability (2–8 years) is low and the reverse for other regions. In addition, statistically significant areas in continuous wavelet spectra show that those timescales of variability are nonstationary (Fig. 4c), with temporal changes in terms of amplitude discriminating the different regions. For instance, southwestern regions are characterized by quasi-continuous significant annual variability until the late 1980s, while other watersheds show sparsely significant annual variability (Fig. 4c). Similarly, although there is significant interannual variability in all watersheds from the late 1980s, during this period, southwestern and southeastern regions do not show significant interannual variations (Fig. 4c). After removing the ≤1-year timescales (i.e., the seasonal cycle) focusing on interannual timescales, significant variability at 2- (southwestern) and 5-year (southwestern and southeastern) timescales emerge for the southern regions. The largest variations (i.e., colored circles) occur over shorter periods of time than in other regions (Fig. 5b).

Figure 5Interannual precipitation time frequency variability in France. (a) Global wavelet spectra of homogeneous regions. (b) Wavelet spectra of homogeneous regions. For clarity, only timescales and time locations that are 95 % statistically significant and with the largest variability are shown (colored circles).

In summary, different regions with coherent precipitation variability are identified and are characterized by three timescales of variability, i.e., intraseasonal, annual, and interannual. The amplitude of those timescales of variability, however, differs in time and over the French territory. Mediterranean regions (southwestern and southeastern) have comparatively weaker interannual variability as compared to annual timescales. The differences between regions are both dependent on the local expression of the climate forcing and watershed characteristics. Because those physical processes are interacting, studying cross-scale interactions in precipitation brings more insight on the dynamics behind the spectral characteristics of each region (Boé2013; Materia et al.2021; Bellucci et al.2015; Ardilouze et al.2020).

Figure 6Precipitation cross-scale interactions (95 % significance level). The driving timescale is on the horizontal axis and the driven on the vertical axis (i.e., the timescale where the x phase has a causal relationship with the phase/amplitude of the driven timescale y). The lower (upper) half of the graph, below (above) the diagonal, shows timescales acting on smaller (larger) timescales. (a) Phase–phase causality. (b) Phase–amplitude causality.

3.1.2 Cross-scale interactions

Figure 6 shows cross-scale interactions for each cluster of precipitation variability (see Fig. 4).

Northeastern, southeastern, north-central, northwestern, and central-eastern regions all show the phase of a 5–8-year variability driving the variability of smaller timescales (Fig. 6a; blue, dark green, red, green, and black; lower half of the graph). This cross-scale interaction is, however, more pronounced in the northeastern and southeastern regions (Fig. 6a). Similarly, eastern regions exclusively show 5–8 to 2–4-year interactions, while other regions show the self-interacting 5–8-year variability (Fig. 6a). The upper half of the graph, which refers to the higher frequency driving the lower-frequency variability, is populated by north-central, southeastern, northwestern, and northeastern regions (Fig. 6a; red, dark green, green and blue). The southeastern region shows cascade phase–phase interactions, (i.e., from 2–3 and 5–4 to 6–5 years; Fig. 6a, dark green). In addition, both southeastern and northwestern regions show mirror interactions with their lower-half counterparts, e.g., 5–6 to 4–5 years (Fig. 6a; dark green and green mirror patches about the diagonal). We also note that phase–phase interactions are very weak over the southwestern regions and absent in the central-western regions.

Phase–amplitude interactions are presented in Fig. 6b. The lower half of the graph, which refers to the lower frequency driving the higher-frequency variability, shows 5–8- to 2–4-year interactions for the western and north-central regions (Fig. 6b; pink, yellow, green, and red). The central-eastern regions are also showing the lower-frequency variability driving the higher-frequency variability but between 8- and 6-year variability (Fig. 6b). Notably, the northwestern region is the only one with cross-scale interactions driving the annual cycle (Fig. 6a; green). In the upper half of the graph, which refers to higher frequency driving lower-frequency variability, we only find north-central and northeastern regions showing 2–4- to 4- and 3–4- to 7–8-year phase–amplitude interactions (Fig. 6b; blue and red). Note that north-central, northeastern, and central-eastern regions show phase–amplitude and phase–phase interactions at very similar timescales (Fig. 6a and b; red, blue, and black), while timescales of phase–amplitude and phase–phase interactions do not match in central-western, northwestern, and southwestern regions (Fig. 6a and b; pink, green, and yellow). Regions to the east, thus, appear to have both phase–phase and phase–amplitude interactions at the same timescales, while western regions are more characterized by phase–amplitude interactions.

The precipitation cross-scale interactions can be of different forms, namely phase–phase, phase–amplitude, uni- or bi-directional, and from lower to higher timescales and vice versa. The presence of cross-scale interactions seems to be tied to specific spatial locations, suggesting different internal dynamics over the different regions of homogeneous precipitation variability. Interestingly, cross-scale interactions tend to converge toward specific timescales, notably 2–4 and 5–8 years, which were linked to ocean–atmosphere variability, such as the North Atlantic Oscillation, in previous hydroclimate studies over France (Feliks et al.2011; Fritier et al.2012; Dieppois et al.2016; Massei et al.2017). In addition, the presence of mirror interactions also indicate strong bidirectional negative feedback.

Figure 7Clustering of temperature time frequency variability in France. (a) Classification map of the watersheds. Pie chart slices show the three highest probability memberships. (b) Global wavelet spectra of regions. (c) Wavelet spectra of homogeneous regions. For clarity, only timescales and time locations that are 95 % statistically significant and with the largest variability are shown (colored circles).

3.2 Temperature

3.2.1 Time frequency patterns

In temperature, the following nine regions with homogeneous time frequency patterns are identified (Fig. 7a): northwestern high (pink), northwestern low (black), northeastern (blue), central eastern (red), central western (green), southeastern high (yellow), southeastern low (brown), southwestern high (dark green), and southwestern low (purple). Fuzzy clustering shows that watersheds typically converge toward singular clusters, defining highly coherent regions (Fig. 7a). This is, however, not true for the central-western region, which is characterized by a mix of spectral characteristic defining other regions (see the red, green, black, yellow, and purple pie charts; Fig. 7a).

Using monthly data, temperature is primarily varying on an annual timescale, with very similar amplitudes for all regions (Fig. 7b and c). Since the dominant annual variability masks the other timescales, we use the annual time step to study interannual variability (Fig. 8a and b). Focusing on this interannual variability, significant temperature variations indeed emerge at 2–4- and 5–8-year timescales and show different timings and amplitudes over the different regions (Fig. 8a and b). All regions show 5–8-year variability, but, compared to northern regions, southern regions show significantly stronger variations on 2–4-year timescale (Fig. 8a). Similarly, while stronger 2–4-year variability occurs in the 1980s and 1990s in the southwestern low region, other regions show significant 2–4-year variability from the 2000s only (brown; Fig. 8b).

Figure 8Interannual temperature time frequency variability in France. (a) Global wavelet spectra of homogeneous regions. (b) Wavelet spectra of homogeneous regions. For clarity, only timescales and time locations that are 95 % statistically significant and with the largest variability are shown (colored circles).

Figure 9Temperature cross-scale interactions (95 % significance level). The driving timescale is on the horizontal axis and the driven on the vertical axis (i.e., the timescale x phase has a causal relationship with the phase/amplitude of the driven timescale y). Lower (upper) half of the graph, below (above) the diagonal, shows timescales acting on smaller (larger) timescales. (a) Phase–phase causality. (b) Phase–amplitude causality.

3.2.2 Cross-scale interactions

Figure 9 shows cross-scale interactions for each cluster of temperature variability identified in Fig. 8a.

For temperature, phase–phase interactions are mostly concentrated in the upper half of the graph, which refers to higher frequency modulating lower frequency (Fig. 9a). Notably, a 2–6- to 6–8-year phase–phase interaction appears more pronounced over northern regions (Fig. 9a; blue, red, pink, and black). The central-western region shows similar phase–phase interactions, but at 1–3- to 4–6-year timescales (Fig. 9a; green). In the lower half of the graph, which refers to lower frequency modulating higher frequency, interactions are found at very similar timescales, but at slightly higher frequency, for all regions (e.g., 2–5- to 1–4-year variability; Fig. 9a). Temperature in the southwestern low region, however, show slightly different characteristics with phase–phase interactions between lower and higher frequency occurring between 7–8- and 3–4- and 7–8- to 3–4-year variability (Fig. 9a; purple).

Temperature phase–amplitude interactions are mostly acting on the 3–4-year timescale for all regions (Fig. 9b). In particular, in temperature, more pronounced phase–amplitude interactions are found over the southwestern low region (Fig. 9b; purple), which is consistent with previous studies on phase–amplitude interactions in European temperature (Palus2014; Jajcay et al.2016). Over southwestern regions, temperature, however, shows both 3–8- and 3–4- and 2–4- and 4–7-year phase–amplitude interactions (Fig. 9b; brown and purple). Furthermore, it should be noted that temperature variability interactions occur between very similar timescales over a number of regions (Fig. 9b; pink, red, yellow, and purple). According to Palus (2014), interactions between very similar timescales, or the same timescales, can only occur if, at least, two processes are present.

As for precipitation, in temperature, phase–phase and phase–amplitude cross-scale interactions are region dependent and can be uni- or bi-lateral. However, in temperature, most phase–phase interactions occur from higher- to lower-frequency variability, while phase–amplitude interactions tend to occur from lower- to higher-frequency variability. Similarly, while timescales of variability that are involved for phase–phase and phase–amplitude interactions are very similar in precipitation, they differ largely in temperature (Fig. 9b). This suggests that, in temperature, the processes driving phase–phase and phase–amplitude cross-scale interactions are different. It also suggests that the processes driving cross-scale interactions are different in temperature and in precipitation.

Figure 10Clustering of discharge time frequency variability in France. (a) Classification map of the watersheds. Pie chart slices show the three highest probability memberships. (b) Global wavelet spectra of homogeneous regions. (c) Wavelet spectra of homogeneous regions. For clarity, only timescales and time locations that are 95 % statistically significant and with the largest variability are shown (colored circles).

3.3 Discharge

3.3.1 Time frequency patterns

The following six regions with homogeneous time frequency patterns are identified in discharge (Fig. 10a): northwestern (black), northeastern (blue), north central (red), central western (green), southeastern (yellow) and southwestern (pink). However, several watersheds, especially in the south, show memberships in multiple regions, suggesting lower spatial coherence in discharge than in precipitation and temperature. Lower spatial coherence, however, could mostly be explained by (i) mixing of solid and liquid precipitation in driving discharge variability in the Alps and (ii) the local heterogeneity of precipitation due to convective dynamics in the Pyrenees (Gottardi et al.2008; Büntgen et al.2008; Hermida et al.2015). Nevertheless, the number of significant homogeneous regions is lower in discharge than in precipitation and temperature, and northern regions are particularly coherent.

Using monthly data, discharge is mainly varying on annual timescales, as determined through the global wavelet spectra (Fig. 10b). In addition, unlike other regions, southeastern watersheds show significant intraseasonal variability (Fig. 10b). Continuous wavelet spectra show that both annual and intraseasonal variability can be nonstationary, with temporal changes in terms of amplitude discriminating the different homogeneous regions of discharge variability (Fig. 10c). For instance, annual variability is only significant for specific periods in the southeastern watersheds, while other regions show quasi-continuous significant annual variability (Fig. 10c). Similarly, in the southeastern region, intraseasonal discharge variability sparsely appears significant from the 1980s, while they are absent in other regions (Fig. 10b).

Figure 11Interannual discharge time frequency variability in France. (a) Global wavelet spectra of homogeneous regions. (b) Wavelet spectra of homogeneous regions. For clarity, only timescales and time locations that are 95 % statistically significant and with the largest variability are shown (colored circles).

After removing the seasonality, and focusing on interannual variability, northeastern watersheds stand out as having continuous significant interannual variability throughout the time series, with 4–5- and 5–8-year variability before and after the 1990s, respectively (Fig. 11b). Southeastern and southwestern regions also stand out, as they show 2–4-year variability in the mid-1970s and 2000s (Fig. 11b; yellow and pink). In addition, southeastern regions do not show significant variability in discharge at timescales greater than 4 years (Fig. 11a and b).

Different coherent regions are thus identified for discharge variability. In addition, these homogeneous regions correspond well with regions identified in precipitation and temperature variability. As in precipitation and temperature, those regions seem strongly impacted constrained temperature, and southern regions, which may appear more complex in term of climate and its link to land surface processes, appear much less spatially coherent in discharge.

Figure 12Discharge cross-scale interactions (95 % significance level). The driving timescale is on the horizontal axis and the driven on the vertical axis (i.e., the timescale x phase has a causal relationship with the phase/amplitude of the driven timescale y). Lower (upper) half of the graph, below (above) the diagonal, shows timescales acting on smaller (larger) timescales. (a) Phase–phase causality. (b) Phase–amplitude causality.

3.3.2 Cross-scale interactions

An important question concerning discharge cross-scale interactions is whether interactions found in either precipitation or temperature are also present in discharge. Phase–phase interactions that were found in precipitation are also identified in discharge, in particular over the northeastern, southeastern, and northwestern regions (blue, yellow, and black; Figs. 6a and 12a). Phase–phase interactions that were identified in temperature are much less evident (Figs. 9a and 12a). It should also be noted that many small patches, describing phase–phase interactions in precipitation and temperature, are systematically not transferred to discharge variability (Figs. 6a, 9a, 12a). Instead, discharge variability seems to exclusively preserve large patches of phase–phase interactions (Figs. 6a, 9a, 12a), suggesting that catchment properties are modulating the climatic signals (i.e., precipitation and temperature). Such filtering of climate signals is even more pronounced in certain regions, such as the north central, where phase–phase interactions are absent in discharge (Fig. 12a) but were identified in precipitation and temperature (Figs. 6a and 9a).

More importantly, there is no phase–amplitude interaction in discharge (Fig. 12b). This points out that watershed properties modulate the interacting processes in precipitation and temperature. Because our data set is mostly composed of low groundwater support, those modulations are unlikely to result from the water table, especially as phase–phase interactions are inherited from precipitation. In addition, further analysis at the Paris Austerlitz gauging station, which includes very large groundwater support, reveals the same absence of phase–amplitude interaction in discharge (not shown; Flipo et al.2020). Possible explanations include the frequency partitioning of watershed compartments or integration process along the river network breaks any spatial connection and thus smooths out and flattens phase–amplitude interactions (Schuite et al.2019).

Cross-scale interactions are only of phase–phase nature in discharge. All phase–phase interactions in discharge seem to be primarily related to precipitation, even though the strong correlations between rainfall and temperature makes it difficult to detect. However, differences between regions of homogeneous discharge variability are very similar to those detected in precipitation. Further work is, however, needed to understand why phase–amplitude cross-interactions are absent in discharge variability. Catchment properties appear to involve positive rather than negative feedback, thus resulting in a loss in phase–amplitude interactions.

4 Discussions and conclusions

4.1 Spatial variability of homogeneous hydroclimate variability in France

As recommended by Blöschl et al. (2019), characterizing the different scales of spatial and temporal variability, as well as their interactions, remains one of the most important challenges in hydrology. In this study, we unraveled homogeneous regions of hydroclimate variability in France, accounting for nonstationarity and nonlinearity, bringing additional information over previous, regime-based, classifications in France or elsewhere (Champeaux and Tamburini1996; Bower and Hannah2002; Sauquet et al.2008; Snelder et al.2009; Joly et al.2010; Gudmundsson et al.2011). This was achieved through a clustering analysis based on time frequency patterns of precipitation, temperature, and discharge variability over 152 watersheds. We then studied the spatiotemporal characteristics of each homogeneous region, including characteristic timescales of hydroclimate variability (i.e., precipitation, temperature, and discharge) and cross-scale interactions.

Our study reveals different coherent regions of precipitation, temperature and discharge variability. Yet, some watersheds are characterized by a mix of spectral characteristics from surrounding regions or regions with the same topographic characteristics. Those coherent regions are homogeneously distributed over France in precipitation and discharge but show larger discrepancies in term of spatial extension in temperature. According to previous clustering of hydroclimate variability over France, northern regions are more homogeneous than what was found here (Champeaux and Tamburini1996; Sauquet et al.2008; Snelder et al.2009; Joly et al.2010) and show lower spatial coherence. In particular, here, we demonstrate that both the amplitude and timings of the different timescales of hydroclimate variability differentiate the regions, highlighting the need for accounting for nonstationary behavior in global to regional hydroclimate study. Overall, hydroclimate variability displays intraseasonal (<1 year), annual (1 year), and interannual (2–4 and 5–8 years) timescales. Our results, which were focused on the French territory, are therefore consistent with timescales of variability identified over the world's major rivers (Labat2006).

The timescales identified in this study have been shown to be important in climate processes, such as the North Atlantic Oscillation or the Gulf Stream front (Massei et al.2007; Feliks et al.2010; O'Reilly et al.2017). Their interactions with watershed characteristics likely leads to their local expression with local processes, which play an important role in feedback mechanisms dampening or enhancing how the climate variability is expressed at the local scale (Haslinger et al.2021; Materia et al.2021; Bellucci et al.2015).

4.2 Cross-scale interactions

Feedback mechanisms can occur between any physical processes of the hydroclimate system, and identifying or attributing the nature of these processes is an intractable issue using observational data. Nevertheless, we can use the mandatory conditions for cross-scale interactions to discuss the processes that are potentially at play (Fig. 3). In precipitation, cross-scale interactions involve lower-frequency timescales driving higher-frequency timescales. North Atlantic climate variability encompasses various processes, such as North Atlantic Oscillation or sea surface temperature anomalies, that drive climate variability (Feliks et al.2010; O'Reilly et al.2016). Thus, moisture advection from the North Atlantic area could potentially act as a positive feedback forcing. Moisture advection has indeed been shown to impact western Europe precipitation, especially in wintertime (Sun et al.2020; O'Reilly et al.2017). Zonal moisture advection is only forcing precipitation variability when the region is not affected by blocking weather regimes (Haslinger et al.2019, 2021). Furthermore, vegetation, temperature and soil moisture, which are themselves interacting with each other, can act as a dampening forcing, thus dampening the precipitation. The precipitation–temperature, precipitation–soil moisture and precipitation–vegetation feedback have indeed been shown to reach a negative sign, depending on prior state of the soil (Liu et al.2006; Berg et al.2015; Liu et al.2006). However, the sign of temperature–soil moisture–vegetation feedback on precipitation has been shown to be spatially dependent at the global scale. For instance, while temperature and soil moisture have large effects in western Europe, vegetation feedback is stronger and mostly of a positive sign in northern Europe (Woodward et al.1998; Liu et al.2006; Yang et al.2018). In our results, the southeastern region shows interannual phase–phase interactions (Fig. 6a), which in contradiction with recent literature. For instance, in the Mediterranean region, Ardilouze et al. (2020) found no negative soil moisture precipitation feedback for interannual scales; however, the authors used two climate models to simulate soil moisture sensitivity to precipitation forcing and note that this variability is much larger in century-long reanalysis, such as NOAA's 20CR. For other regions, interannual negative soil moisture feedback was found by Boé (2013), while Sejas et al. (2014) found negative ocean–land temperature differences in precipitation feedback. Similar results were found in Bellucci et al. (2015), where interactions between compartments of the atmospheric circulation at intraseasonal timescale were found to produce significant interannual variations.

Interannual temperature variability is tied to both the soil state and atmospheric circulation, but that relation is location dependent. Large-scale patterns, such as the North Atlantic Oscillation, are shown to be source of both interannual precipitation and temperature variability, especially during wintertime, including for southwestern France (Pepin and Kidd2006; O'Reilly et al.2016). At more local scales, sea surface temperature anomalies have been shown to interact with near-surface air temperature through sea–land heat exchanges regions (Lambert et al.2011; Sejas et al.2014; Zveryaev2015). Soil moisture and evapotranspiration demand can enhance or dampen near-surface temperature variability (Miralles et al.2012; Materia et al.2021). Here, in temperature, phase–phase interactions are particularly interesting because they arise from higher-frequency timescales driving lower-frequency timescales (Fig. 9a). As shown by Peters et al. (2004, 2007), higher-frequency processes can spread to lower-frequency ones by the means of intermediate timescale processes. High-frequency soil moisture enhancing lower-frequency large-scale circulation may explain temperature cross-scale interactions.

Regarding cross-scale interactions in discharge variability, the absence of phase–amplitude was particularly interesting. As our discharge data set is mostly composed of low groundwater support, the absence of phase–amplitude interactions is unlikely to result from the water table, especially as phase–phase interactions are inherited from precipitation. To test this hypothesis, we computed cross-scale interactions on the gauging station at Paris Austerlitz, which was not included in our original data set, as it shows large groundwater support and anthropogenic influence. Results at Paris Austerlitz are consistent with other regions and do not show any phase–amplitude interactions (not shown). As it has been shown that spatial heterogeneity (in the variable dynamics) favors cross-scale interactions, one possible explanation is that converging of runoff into the main drain cancels that spatial heterogeneity and, thus, phase–amplitude variability (Peters et al.2007).

In this study, we interpreted cross-scale interactions based on the mandatory structure for such interactions to arise, identified interacting timescales, and compared cross-scale interactions in both precipitation, temperature, and discharge. Dedicated studies are needed to explore, in depth, the drivers of those interactions, as feedback mechanisms are complex, and likely different for each variable, even though phase–phase interactions in discharge clearly show the signature of those identified in precipitation.

4.3 Conclusion

Those findings allow for a better identification of climate-deterministic processes controlling hydroclimate variability, notably using cross-scale analysis, which could help in identifying more robust climate drivers. For instance, it is important to discriminate pure climate influence from climate–land processes interactions. This has large implications for seamless hydrological predictions based on climate information, as only some parts of the climate signals are transferred to discharge systems. Thus, causal cross-scale relationships could be used to inform and improve existing seasonal to multi-year seamless forecasting for hydrological variability, including extremes (e.g., flood and drought). Preliminary work in this direction was recently proposed by Jajcay et al. (2018), who developed a composite binning method enabling one to forecast a particular time series based on conditional phase of another. Similarly, it would be of crucial importance to determine whether hydrological models, which are commonly used in climate-impact assessments, are reproducing the filtering processes induced by the catchment properties and then identify those (Ducharne et al.2020). Long-term hydroclimate variability only represents a fraction of the total variability; however, strong interactions between high- and low-frequency variations have been highlighted. Those interactions are both spatial and temporal (Feliks et al.2016). Owing to the recent addition of long-term, high spatial resolution hydroclimate data sets (e.g., Fyre reconstructions; Devers et al.2020, 2021), it is now possible to apply the clustering and cross-scale analyses to better characterize the effects that long-term hydroclimate variability (e.g., multi-decadal) has on smaller timescales. The methodology presented in this work can enable deeper analyses than those based on correlations, which may overlook some important hydroclimate processes.

Code availability

The code used for this study is available at (Fossa2021). The pyclits Python code used for cross-scale interactions is available at (last access: 17 March 2021; Jajcay and Lavicka2018).

Data availability

The Safran precipitation and temperature data set can be obtained from (last access: 1 October 2015). Discharge data are available from (last access: 15 October 2015 Ministère de la Transition Energetique2015).


The figures presented in this article are available in a colorblind-friendly palette in the Supplement. The supplement related to this article is available online at:

Author contributions

MFos, NM, and BD conceptualized the study. MFos took responsibility for the methodology, software, formal analysis, investigation, original draft preparation, and visualization. MFos, BD, NM, and MFou validated the study. JPV collected the resources, and MFos and JPV curated the data. MFos, BD, NM, and JPV reviewed and edited the paper.

Competing interests

The authors declare that they have no conflict of interest.


Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.


The authors would like to thank Nikola Jajcay, for his support with the adaptation of the Conditional Mutual Information algorithm (pyclits) to our context. The authors would also like to thank the Centre Régional Informatique et d'Applications Numériques de Normandie (CRIANN), for providing the HPC environment needed for all computations. This research work is part of a contribution to the EURO-Friend group 2.

Review statement

This paper was edited by Miriam Coenders-Gerrits and reviewed by two anonymous referees.


Anctil, F. and Coulibaly, P.: Wavelet Analysis of the Interannual Variability in Southern Québec Streamflow, J. Climate, 17, 163–173, 2004. a

Ardilouze, C., Materia, S., Batté, L., Benassi, M., and Prodhomme, C.: Precipitation response to extreme soil moisture conditions over the Mediterranean, Clim. Dynam.,, 2020. a, b, c

Bellucci, A., Haarsma, R., Bellouin, N., Booth, B., Cagnazzo, C., van den Hurk, B., Keenlyside, N., Koenigk, T., Massonnet, F., Materia, S., and Weiss, M.: Advancements in decadal climate predictability: The role of nonoceanic drivers, Rev. Geophys., 53, 165–202, 2015. a, b, c, d

Berg, A., Lintner, B. R., Findell, K., Seneviratne, S. I., van den Hurk, B., Ducharne, A., Chéruy, F., Hagemann, S., Lawrence, D. M., Malyshev, S., Meier, A., and Gentine, P.: Interannual Coupling between Summertime Surface Temperature and Precipitation over Land: Processes and Implications for Climate Change, J. Climate, 28, 1308–1328, 2015. a

Blöschl, G., Bierkens, M. F. P., Chambel, A., Cudennec, C., Destouni, G., Fiori, A., Kirchner, J. W., McDonnell, J. J., Savenije, H. H. G., Sivapalan, M., Stumpp, C., Toth, E., Volpi, E., Carr, G., Lupton, C., Salinas, J., Széles, B., Viglione, A., Aksoy, H., Allen, S. T., Amin, A., Andréassian, V., Arheimer, B., Aryal, S. K., Baker, V., Bardsley, E., Barendrecht, M. H., Bartosova, A., Batelaan, O., Berghuijs, W. R., Beven, K., Blume, T., Bogaard, T., Borges de Amorim, P., Böttcher, M. E., Boulet, G., Breinl, K., Brilly, M., Brocca, L., Buytaert, W., Castellarin, A., Castelletti, A., Chen, X., Chen, Y., Chen, Y., Chifflard, P., Claps, P., Clark, M. P., Collins, A. L., Croke, B., Dathe, A., David, P. C., de Barros, F. P. J., de Rooij, G., Di Baldassarre, G., Driscoll, J. M., Duethmann, D., Dwivedi, R., Eris, E., Farmer, W. H., Feiccabrino, J., Ferguson, G., Ferrari, E., Ferraris, S., Fersch, B., Finger, D., Foglia, L., Fowler, K., Gartsman, B., Gascoin, S., Gaume, E., Gelfan, A., Geris, J., Gharari, S., Gleeson, T., Glendell, M., Gonzalez Bevacqua, A., González-Dugo, M. P., Grimaldi, S., Gupta, A. B., Guse, B., Han, D., Hannah, D., Harpold, A., Haun, S., Heal, K., Helfricht, K., Herrnegger, M., Hipsey, M., HlaváĊiková, H., Hohmann, C., Holko, L., Hopkinson, C., Hrachowitz, M., Illangasekare, T. H., Inam, A., Innocente, C., Istanbulluoglu, E., Jarihani, B., Kalantari, Z., Kalvans, A., Khanal, S., Khatami, S., Kiesel, J., Kirkby, M., Knoben, W., Kochanek, K., Kohnová, S., Kolechkina, A., Krause, S., Kreamer, D., Kreibich, H., Kunstmann, H., Lange, H., Liberato, M. L. R., Lindquist, E., Link, T., Liu, J., Loucks, D. P., Luce, C., Mahé, G., Makarieva, O., Malard, J., Mashtayeva, S., Maskey, S., Mas-Pla, J., Mavrova-Guirguinova, M., Mazzoleni, M., Mernild, S., Misstear, B. D., Montanari, A., Müller-Thomy, H., Nabizadeh, A., Nardi, F., Neale, C., Nesterova, N., Nurtaev, B., Odongo, V. O., Panda, S., Pande, S., Pang, Z., Papacharalampous, G., Perrin, C., Pfister, L., Pimentel, R., Polo, M. J., Post, D., Prieto Sierra, C., Ramos, M.-H., Renner, M., Reynolds, J. E., Ridolfi, E., Rigon, R., Riva, M., Robertson, D. E., Rosso, R., Roy, T., Sá, J. H. M., Salvadori, G., Sandells, M., Schaefli, B., Schumann, A., Scolobig, A., Seibert, J., Servat, E., Shafiei, M., Sharma, A., Sidibe, M., Sidle, R. C., Skaugen, T., Smith, H., Spiessl, S. M., Stein, L., Steinsland, I., Strasser, U., Su, B., Szolgay, J., Tarboton, D., Tauro, F., Thirel, G., Tian, F., Tong, R., Tussupova, K., Tyralis, H., Uijlenhoet, R., van Beek, R., van der Ent, R. J., van der Ploeg, M., Van Loon, A. F., van Meerveld, I., van Nooijen, R., van Oel, P. R., Vidal, J.-P., von Freyberg, J., Vorogushyn, S., Wachniew, P., Wade, A. J., Ward, P., Westerberg, I. K., White, C., Wood, E. F., Woods, R., Xu, Z., Yilmaz, K. K. and Zhang, Y.: Twenty-three Unsolved Problems in Hydrology (UPH) – a community perspective, Hydrol. Sci. J., 0, 1–33,, 2019. a, b, c

Boé, J.: Modulation of soil moisture–precipitation interactions over France by large scale circulation, Clim. Dynam., 40, 875–892, 2013. a, b

Bower, D. and Hannah, D. M.: Spatial and temporal variability of UK river flow regimes, IAHS-AISH Publication, pp. 457–466, 2002. a

Brigode, P., Génot, B., Lobligeois, F., and Delalgue, O.: Summary sheets of watershed-scale hydroclimatic observed data for France, sheets of watershed-scale hydroclimatic observed data for France, available at: (last access: 2 July 2021), 2020. a

Büntgen, U., Frank, D., Grudd, H., and Esper, J.: Long-term summer temperature variations in the Pyrenees, Clim. Dynam., 31, 615–631, 2008. a

Caillouet, L., Vidal, J.-P., Sauquet, E., Devers, A., and Graff, B.: Ensemble reconstruction of spatio-temporal extreme low-flow events in France since 1871, Hydrol. Earth Syst. Sci., 21, 2923–2951,, 2017. a

Champeaux, J.-L. and Tamburini, A.: Zonage climatique de la France à partir des séries de précipitations (1971–1990) du réseau climatologique d’Etat, La météorologie, 8, 44–54, 1996. a, b

Christophe Bouton, P. H. (Ed.): Time of Nature and the Nature of Time: Philosophical Perspectives of Time in Natural Sciences, vol. 326 of (Boston Studies in the Philosophy and History of Science, Springer, 2017. a

Coulibaly, P. and Burn, D. H.: Wavelet analysis of variability in annual Canadian streamflows, Water Resour. Res., 40, 1–14,, 2004. a, b, c

Devers, A., Vidal, J. P., Lauvernet, C., Graff, B., and Vannier, O.: A framework for high-resolution meteorological surface reanalysis through offline data assimilation in an ensemble of downscaled reconstructions, Q. J. Roy. Meteor. Soc., 146, 153–173,, 2020. a

Devers, A., Vidal, J.-P., Lauvernet, C., and Vannier, O.: FYRE Climate: a high-resolution reanalysis of daily precipitation and temperature in France from 1871 to 2012, Clim. Past, 17, 1857–1879,, 2021. a

Dieppois, B., Durand, A., Fournier, M., and Massei, N.: Links between multidecadal and interdecadal climatic oscillations in the North Atlantic and regional climate variability of northern France and England since the 17th century, J. Geophys. Res.-Atmos., 118, 4359–4372,, 2013. a, b

Dieppois, B., Lawler, D. M., Slonosky, V., Massei, N., Bigot, S., and Fournier, M.: Multidecadal climate variability over northern France during the past 500 years and its relation to large-scale atmospheric circulation, Int. J. Climatol., 36, 4679–4696,, 2016. a, b, c

Ducharne, A., Arboleda-obando, P., and Cheruy, F.: Effets de l’humectation des sols par les nappes sur la trajectoire du changement climatique dans le bassin de la Seine et en Europe, Tech. rep., PIREN-Seine, Paris, 2020. a

Dunn, J. C.: A Fuzzy Relative of the ISODATA Process and Its Use in Detecting Compact Well-Separated Clusters A Fuzzy Relative of the ISODATA Process and Its Use in Detecting Compact Well-Separated Clusters, J. Cybernet., 3, 37–57, 1973. a

Ebisuzaki, W.: A Method to Estimate the Statistical Significance of a Correlation When the Data Are Serially Correlated, J. Climate, 10, 2147–2153, 1997. a

Feliks, Y., Ghil, M., and Robertson, A. W.: Oscillatory Climate Modes in the Eastern Mediterranean and Their Synchronization with the North Atlantic Oscillation, J. Climate, 23, 4060–4079, 2010. a, b

Feliks, Y., Ghil, M., and Robertson, A. W.: The atmospheric circulation over the North Atlantic as induced by the SST field, J. Climate, 24, 522–542,, 2011. a

Feliks, Y., Robertson, A. W., and Ghil, M.: Interannual variability in north Atlantic weather: Data analysis and a quasi565 geostrophic model, J. Atmos. Sci., 73, 3227–3248,, 2016. a

Flipo, N., Gallois, N., Labarthe, B., Baratelli, F., Viennot, P., Schuite, J., Rivière, A., Bonnet, R., and Boé, J.: Pluri-annual Water Budget on the Seine Basin: Past, Current and Future Trends, in: Handbook of environmental chemistry, Springer,, 2020. a

Fossa, M.: ManuelFossa/Hess-2021-81: HESS-2021-81 (CMI), Zenodo [code],, 2021. a

Fritier, N., Massei, N., Laignel, B., Durand, A., Dieppois, B., and Deloffre, J.: Links between NAO fluctuations and inter-annual variability of winter-months precipitation in the Seine River watershed (north-western France), Comptes Rendus – Geoscience, 344, 396–405,, 2012. a

Ge, Z.: Significance tests for the wavelet power and the wavelet power spectrum, Ann. Geophys., 25, 2259–2269,, 2007. a, b, c

Gentine, P., Troy, T. J., Lintner, B. R., and Findell, K. L.: Scaling in Surface Hydrology: Progress and Challenges, J. Contemp. Water Res. Educ., 147, 28–40,, 2012. a, b

Giuntoli, I., Renard, B., Vidal, J.-P., and Bard, A.: Low flows in France and their relationship to large-scale climate indices, J. Hydrol., 482, 105–118,, 2013. a, b, c

Gottardi, F., Obled, C., Gailhard, J., and Paquet, E.: Régionalisation des précipitations sur les massifs montagneux français à l'aide de régressions locales et par types de temps, Climatologie, 5, 7–25, 2008. a

Granger, C. W. J.: Investigating Causal Relations by Econometric Models and Cross-spectral Methods, Econometrica, 37, 424–438, 1969. a

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. a

Gudmundsson, L., Tallaksen, L. M., and Stahl, K.: Spatial cross-correlation patterns of European low, mean and high flows, Hydrol. Process., 25, 1034–1045,, 2011. a, b, c

Hannachi, A., Straus, D. M., Franzke, C. L. E., Corti, S., and Woollings, T.: Low-frequency nonlinearity and regime behavior in the Northern Hemisphere extratropical atmosphere: NONLINEARITY AND REGIME BEHAVIOR, Rev. Geophys., 55, 199–234, 2017. a

Hannaford, J., Lloyd-Hughes, B., Prudhomme, C., Parry, S., Keef, C., and Rees, G.: The Spatial Coherence of European Droughts – Final Report protecting and improving the environment in England and Wales, Tech. rep., Environment Agency's Science Programme, Environment Agency, Rio House, Waterside Drive, Aztec West, Almondsbury, Bristol, BS32 4UD, 2009. a

Haslinger, K., Hofstätter, M., Kroisleitner, C., Schöner, W., Laaha, G., Holawe, F., and Blöschl, G.: Disentangling Drivers of Meteorological Droughts in the European Greater Alpine Region During the Last Two Centuries, J. Geophys. Res.-Atmos., 124, 12404–12425, 2019. a

Haslinger, K., Hofstätter, M., Schöner, W., and Blöschl, G.: Changing summer precipitation variability in the Alpine region: on the role of scale dependent atmospheric drivers, Clim. Dynam., 57, 1009–1021, 2021. a, b

Hermida, L., López, L., Merino, A., Berthet, C., García-Ortega, E., Sánchez, J. L., and Dessens, J.: Hailfall in southwest France: Relationship with precipitation, trends and wavelet analysis, Atmos. Res., 156, 174–188, 2015. a

Hubert, P.: Les multifractals, un outil pour surmonter les problèmes d'échelle en hydrologie, Hydrol. Sci. J., 46, 897–905,, 2001. a

Hubert, P., Carbonnel, J. P., and Chaouche, A.: Segmentation Des Séries Hydrométéorologiques – Application à des Séries de Précipitations et de Débits de L'Afrique de l'Ouest, J. Hydrol., 110, 349–367,, 1989. a

IPCC: Climate Change 2007: The Physical Science Basis, Cambridge university press,, 2007. a

IPCC: Climate change 2014. Synthesis report, Cambridge university press,, 2014. a

IPCC: Water Cycle changes, in: IPCC AR6 WGI Full Report, Cambridge University Press, in press, 2021. a

Jajcay, N. and Lavicka, H.: pyCliTS, available at: (last access: 17 March 2021), GitHub [code], 2018. a

Jajcay, N., Hlinka, J., Kravtsov, S., Tsonis, A. A., and Paluš, M.: Time scales of the European surface air temperature variability : The role of the 7–8 year cycle, Geophys. Res. Lett., 43, 1–8,, 2016. a

Jajcay, N., Kravtsov, S., Sugihara, G., Tsonis, A. A., and Paluš, M.: Synchronization and causality across time scales in El Niño Southern Oscillation, npj Climate and Atmospheric Science, 1, 33,, 2018. a, b, c, d

Joly, D., Brossard, T., Cardot, H., Cavailhes, J., Hilal, M., and Wavresky, P.: Les types de climats en France , une construction spatiale, Cybergeo: European Journal of Geography, Document 501,, 2010. a, b

Kaufman, L. and Rousseeuw, P. J.: Finding Groups in Data: An Introduction to Cluster Analysis, Wiley interscience, 368 pp., ISBN13 9780471878766, 1990. a, b

Labat, D.: Non-Linéarité et Non-Stationnarité en Hydrologie Karstique, PhD thesis, INP Toulouse, 2000. a

Labat, D.: Oscillations in land surface hydrological cycle, Earth Planet. Sc. Lett., 242, 143–154,, 2006. a, b, c, d, e

Lambert, F. H., Webb, M. J., and Joshi, M. M.: The relationship between land-ocean surface temperature contrast and radiative forcing, J. Climate, 24, 3239–3256,, 2011. a

Lavers, D., Prudhomme, C., and Hannah, D. M.: Large-scale climate, precipitation and British river flows: Identifying hydroclimatological connections and dynamics, J. Hydrol., 395, 242–255,, 2010. a

Liu, D. and Graham, J.: Simple Measures of Individual Cluster-Membership Certainty for Hard Partitional Clustering, American Statistician, 73, 70–79,, 2018. a, b

Liu, Q., Wen, N., and Liu, Z.: An observational study of the impact of the North Pacific SST on the atmosphere: NORTH PACIFIC IMPACT ON ATMOSPHERE, Geophys. Res. Lett., 33, L18611,, 2006. a, b

Massei, N., Durand, A., Deloffre, J., Dupont, J. P., Valdes, D., and Laignel, B.: Investigating possible links between the North Atlantic Oscillation and rainfall variability in Northwestern France over the past 35 years, J. Geophys. Res.-Atmos., 112, 1–10,, 2007. a, b

Massei, N., Dieppois, B., Hannah, D. M., Lavers, D. A., Fossa, M., Laignel, B., and Debret, M.: Multi-time-scale hydroclimate dynamics of a regional watershed and links to large-scale atmospheric circulation: Application to the Seine river catchment, France, J. Hydrol., 546, 262–275,, 2017. a, b, c, d

Materia, S., Ardilouze, C., Prodhomme, C., Donat, M. G., Benassi, M., Doblas-Reyes, F. J., Peano, D., Caron, L.-P., Ruggieri, P., and Gualdi, S.: Summer temperature response to extreme soil water conditions in the Mediterranean transitional climate regime, Clim. Dynam., 2021. a, b, c, d

McGregor, G.: Hydroclimatology, modes of climatic variability and stream flow, lake and groundwater level variability: A progress report, Prog. Phys. Geogr., 41, 496–512,, 2017. a

Ministère de la Transition Energetique: banque hydro, Ministère de la Transition Energetique [data set], available at:, last access: 15 October 2015. a

Miralles, D. G., van den Berg, M. J., Teuling, A. J., and de Jeu, R. A. M.: Soil moisture-temperature coupling: A multiscale observational analysis, Geophys. Res. Lett., 39, L21707,, 2012. a

Monti, S., Tamayo, P., Mesirov, J., and Golub, T.: Consensus clustering: A resampling-based method for class discovery and visualization of gene expression microarray data, Mach. Learn., 52, 91–118,, 2003. a

Moron, V., Robertson, A. W., Ward, M. N., and Camberlin, P.: Spatial coherence of tropical rainfall at the regional scale, J. Climate, 20, 5244–5263,, 2007. a

Nandi, B., Swiatek, P., Kocsis, B., and Ding, M.: Inferring the direction of rhythmic neural transmission via inter-regional phase-amplitude coupling (ir-PAC), Nat. Sci. Rep., 9, 1–13,, 2019. a

Onslow, A. C. E., Jones, M. W., and Bogacz, R.: A canonical circuit for generating phase-amplitude coupling, PLoS ONE, 9, e102591,, 2014. a, b

O'Reilly, C. H., Minobe, S., and Kuwano-Yoshida, A.: The influence of the Gulf Stream on wintertime European blocking, Clim. Dynam., 47, 1545–1567, 2016. a, b

O'Reilly, C. H., Minobe, S., Kuwano‐Yoshida, A., and Woollings, T.: The Gulf Stream influence on wintertime North Atlantic jet variability, Q. J. Roy. Meteor. Soc., 143, 173–183, 2017. a, b

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

Palus, M.: Multiscale atmospheric dynamics: Cross-frequency phase-amplitude coupling in the air temperature, Phys. Rev. Lett., 112, L078702,, 2014. a, b, c, d, e

Pella, H., Lejot, J., Lamouroux, N., and Snelder, T.: The theoretical hydrographical network (RHT) for France and its environmental attributes, Geomorphologie: Relief, Processus, Environnement, 18, 317–336,, 2012. a

Pepin, N. and Kidd, D.: Spatial temperature variation in the Eastern Pyrenees, Weather, 61, 300–310, 2006. a

Peters, D. P. C., Pielke, Sr, 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. a

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. a, b, c

Pikovsky, A., Rosenblum, M., and Kurths, J.: Synchronization. A Universal Concept in Nonlinear Sciences, Cambridge University Press, Cambridge, 2001. a

Rahiz, M. and New, M.: Spatial coherence of meteorological droughts in the UK since 1914, Area, 44, 400–410,, 2012. a

Sauquet, E., Gottschalk, L., and Krasovskaia, I.: Estimating mean monthly runoff at ungauged locations: an application 680 to France, Hydrol. Res., 39, 403–423,, 2008. a, b

Scaife, A. A. and Smith, D.: A signal-to-noise paradox in climate science, npj Climate and Atmospheric Science, 1, 1–8, 2018. a

Schaefli, B., Maraun, D., and Holschneider, M.: What drives high flow events in the Swiss Alps? Recent developments in wavelet spectral analysis and their application to hydrology, Adv. Water Resour., 30, 2511–2525,, 2007. a

Scheffer-Teixeira, R. and Tort, A. B.: On cross-frequency phase-phase coupling between theta and gamma oscillations in the hippocampus, Elife, 5, e20515,, 2016. a

Schuite, J., Flipo, N., Massei, N., Rivière, A., and Baratelli, F.: Improving the Spectral Analysis of Hydrological Signals to Efficiently Constrain Watershed Properties, Water Resour. Res., 55, 4043–4065,, 2019.  a

Sejas, S. A., Albert, O. S., Cai, M., and Deng, Y.: Feedback attribution of the land-sea warming contrast in a global warming simulation of the NCAR CCSM4, Environ. Res. Lett., 9,, 2014. a, b

Şenbabaoǧlu, Y., Michailidis, G., and Li, J. Z.: Critical limitations of consensus clustering in class discovery, Sci. Rep.-UK, 4, 6207,, 2014. a

Sidibe, M., Dieppois, B., Eden, J., Mahé, G., Paturel, J. E., Amoussou, E., Anifowose, B., and Lawler, D.: Interannual to Multi-decadal streamflow variability in West and Central Africa: Interactions with catchment properties and large-scale climate variability, Global Planet. Change, 177, 141–156,, 2019. a

Smith, L., Turcotte, D., and Isacks, B.: Stream flow characterization and feature detection using a discrete wavelet transform, Hydrol. Process., 12, 233–249,<233::AIDHYP573>3.0.CO;2-3, 1998. a

Snelder, T. H., Lamouroux, N., Leathwick, J. R., Pella, H., Sauquet, E., and Shankar, U.: Predictive mapping of the natural flow regimes of France, J. Hydrol., 373, 57–67,, 2009. a, b

Sun, Y., Liu, J.-W., and Xie, S.-P.: North Atlantic Oscillation Effect on Interannual Variability in Winter Precipitation over the Gulf Stream, J. Climate, 33, 6633–6649, 2020. a

Torrence, C. and Compo, G. P.: A practical guide to wavelet analysis, B. Am. Meteorol. Soc., 79, 61–78, 1998. a, b, c, d, e, f

Vidal, J. P., Martin, E., Franchistéguy, L., Baillon, M., and Soubeyroux, J. M.: A 50-year high-resolution atmospheric reanalysis over France with the Safran system, Int. J. Climatol., 30, 1627–1644,, 2010. a, b

Wang, C., Osiński, M., Even, J., and Grillot, F.: Phase-amplitude coupling characteristics in directly modulated quantum dot lasers, Appl. Phys. Lett., 105, L221114,, 2014. a

Wang, L., Zhang, Y., and Feng, J.: On the Euclidean distance of images, IEEE T. Pat. Anal. Mach. Int., 27, 1334–1339,, 2005. a, b, c, d

Woodward, F. I., Lomas, M. R., Betts, R. A., Wagner, A., Mulligan, M., and Hewitt, C. N.: Vegetation-Climate Feedbacks in a Greenhouse World [and Discussion], Philos. Trans. R. Soc. Lond. B Biol. Sci., 353, 29–39, 1998. a

Yang, L., Sun, G., Zhi, L., and Zhao, J.: Negative soil moisture-precipitation feedback in dry and wet regions, Sci. Rep.-UK, 8, 1–9, 2018. a

Zveryaev, I. I.: Seasonal differences in intraseasonal and interannual variability of Mediterranean Sea surface temperature, J. Geophys. Res.-Oceans, 120, 2813–2825, 2015. a

Short summary
Hydro-climate observations (such as precipitation, temperature, and river discharge time series) reveal very complex behavior inherited from complex interactions among the physical processes that drive hydro-climate viability. This study shows how even small perturbations of a physical process can have large consequences on some others. Those interactions vary spatially, thus showing the importance of both temporal and spatial dimensions in better understanding hydro-climate variability.