SPATIO-TEMPORAL AND CROSS-SCALE INTERACTIONS IN HYDROCLIMATE VARIABILITY: A CASE-STUDY IN FRANCE

Understanding how water resources vary at different temporal and spatial scales in response to climate 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 on 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 5 phase-phase and phase-amplitude causal interactions between time scales of each homogeneous region. Three significant time scales of variability are found in precipitation, temperature and discharge: 1 year, 2-4 years and 5-8 years. The magnitude of those time scales 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, we note 10 a sudden decrease in annual temperature variability in the mid 1990s in the Southern half of France. Investigating crossscale 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 timefrequency patterns (occurrence and timing of time scales of variability) were similar between regions, cross-scale interactions are far much complex, differ between regions, and are not systematically transferred from climate variability (precipitation and 15 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-time scale representation of hydro-climate variability in France, and provides unique insight into the complex non-linear dynamics of this variability, and its predictability. 20


Introduction
Hydroclimate variability represents the spatio-temporal evolution of hydrological variables (e.g. discharge, groundwater level) and climate variables (e.g. precipitation and temperature), which are directly impacting hydrological variability. Studying how network of stations identifies near-natural watersheds with long-term high-quality hydrometric data. The period 1968The period -2008 was chosen by Giuntoli et al. (2013) as being the best trade-off in terms of data availability over the different regions. This database was further subset to 152 watersheds by selecting only continuous monthly time series, i.e. without missing values ( Figure 1). Precipitation and temperature data have been estimated from the 8 km-grid Safran surface reanalysis dataset (Vidal et al., 2010), which has been subset to a common period . For this study, precipitation and temperature have been 75 averaged over each watershed area (Caillouet et al., 2017). Figure 2 summarizes the different methods employed in this study.

Continuous wavelet transforms
Time-frequency patterns have first been extracted for each watershed and each variable using continuous wavelet transform (cf. 80 Figure 2, (a)). For any finite energy signal x, it is possible to obtain a time-frequency representation by mapping it to a series of sub-spaces spawned by a generating function, the mother wavelet, and its scaled versions (Torrence and Compo, 1998;Grinsted et al., 2004). The time series is then represented in terms of a given scale and time location. The first subspace is generated by a mother wavelet at scale 1 and its time translations. Then, other subspaces are generated by scaling the mother wavelet up, referring to daughter wavelets, and time translating it. For each scale, one subspace is constructed. Daughter wavelets are 85 usually calculated as: The left hand side (LHS) term is the daughter wavelet of scale a and time translation b at time t. 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 is of the form: LHS term contains the wavelet coefficients, i.e. the coordinates of the signal in each subspace. If the mother wavelet (and hence the daughter wavelets as well) is complex, wavelet coefficients are complex as well. Wavelet coefficients represent the inner product of the signal and daughter wavelet of scale a and time translation b (Centre). The norm of their square is called the wavelet power and represents the amplitude of the oscillation of signal x at scale a and centred on time t. As it is impossible to 95 capture the best resolution in both frequency and time simultaneously, we used here a Morlet mother wavelet (order 6), which offers a good trade-off between detection of scales and localisation of the oscillations in time (Torrence and Compo 1998).
Statistical significance is typically tested against a red noise using Monte Carlo simulations (Torrence and Compo 1998).

Image Euclidean Distance Clustering
We estimated similarities between wavelet spectra of each watersheds and, separately, on each variable. Distances between 100 two-dimensional data, such as maps or wavelet spectra, are estimated using Euclidean distance between pairwise points (pED; i.e computing f 2 (x i , y i )−f 1 (x i , y i )). However, such a procedure has no neighborhood notion, making it impossible to account for globally similar shapes (cf. definition of global and local similarity in 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 two ways (Wang et al., 2005): i) the distance between pixel values is computed not only pairwise, but for 105 all indices; ii) a Gaussian filter, 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 (?).

Fuzzy clustering
Fuzzy clustering has then been used to cluster the different watershed based on their similarities ( Figure 2, (b)). Fuzzy clustering is a soft clustering method (Dunn, 1973). While soft clustering spreads membership over all clusters with varying probability, 110 hard clustering attributes each station one and only one cluster membership. Soft clustering is therefore better-suited when the spatial variability, originating from different stations' characteristics, is smooth, such as in hydroclimate data. 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 New, 2012). As such, several stations tend to show transitional or hybrid patterns, and can potentially be member of different clusters, limiting the robustness 115 of hard clustering procedure (Liu and Graham, 2018).
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 Rousseeuw, 1990). Here, we used the FANNY algorithm (Kaufman and Rousseeuw, 1990), which has been shown to be flexible, and to offer the possibility to adapt the clustering 120 to the data, with optimal performance (Liu and Graham, 2018). In addition, rather than setting the number arbitrarily, in this study, we used an estimation of the optimum number of clusters by first computing a hard clustering method: the consensus clustering (Monti et al., 2003). Thus, the number of clusters providing the best stability (i.e. the minimal changes of 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.

cross-scale interactions
For each variable and each cluster, cross-scale interactions have also been explored for the first time in hydrological sciences ( Figure 2, (c)). Cross-scale interactions refer to phase-phase and phase-amplitude couplings between time scales of a given time series (Scheffer-Teixeira and Tort, 2016;Paluš, 2014). 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 (Granger, 1969;Pikovsky et al., 2001), which 130 refers to information transfer from one part of a signal to another. with each other in a perturbation-dampening (X, Y , respectively), so that f (t) = X − Y (Figure 3a). The interactions between the components occur through the connections C XY and C Y X , with a given strength (here C .. = 2), and this perturbation-135 dampening interaction forms a negative feedback system ( Figure 3a). The connection C XX enables X to grow first before Y dampens it, and both X and Y receive inputs φ X , φ Y from other physical processes (figure 3a). Depending on both the mean and time scales of φ X and φ Y , the strength of C XX , C XY and C Y X , X and Y may show coupled behaviors. For instance, in Figure 3b, every fourth ridges of Y P P (t) is synchronized with a ridge of X(t) (Figure 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 C XX , C XY 140 and C Y X . f P P (t) is the difference between X(t) and Y P P (Figure 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 ( Figure 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 go back to normal when outside of the phase of the driving component.
This is a phase-amplitude interaction. In Figure 3b, Y P A amplitude decreases when X is at its maximum (i.e. when its phase 145 is a ridge, Figure 3c, top and middle panels, respectively). Similarly to the phase-phase interaction, f P A (t) is the difference between X(t) and Y P A (t) (Figure 3c, bottom panel). Because phase amplitude are very dependent on inputs φ X and φ Y , connections between spatially distant physical processes are likely to give rise of phase-amplitude interactions (Nandi et al., 2019).
Following Paluš (2014) and (Jajcay et al., 2018), who compared the most used methods when studying causality, we choose 150 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 It is possible for a single time scale 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, mean and standard deviation as the original time series, as in Ebisuzaki (1997). Paluš (2014) has shown that this number of surrogates is ideal for statistical significance, in the context

Spatio-temporal clustering of hydrological variability
The wavelet transforms corresponding to each set of variables and each watershed have been computed and checked for 165 similarities using IEDC fuzzy clustering to identify and characterize homogeneous regions of hydroclimate variability over France. Cross-scale couplings were then investigated for each homogeneous region.

Time-frequency patterns
Seven regions with homogeneous time-frequency patterns are identified ( Figure 4a): North-western (green), North-eastern 170 (blue), Centre-North (red), Centre-western (pink), Centre-eastern (black), South-western (yellow) and South-eastern (dark green). Based on fuzzy clustering, except for the Northern Alps regions, which show a mix of different regions (i.e. Centreeastern, Centre-western, North-eastern, Centre-North and North-western), most regions are dominated by a spectral characteristic, defining a single member (i.e. cluster), suggesting a greater spatial coherence. Note that Centre-eastern, Centre-western, North-eastern, Centre-North and North-western regions are characterized by mountainous ranges, where there is a significant 175 part of solid precipitation, which may contribute to membership mixing, and weaker spatial coherence in those regions. Similarly, Northern regions are split between oceanic influence (to the west), inland (Centre) and mountainous ranges. Centre and Southern regions are delineated by mountain and valleys spatial extent (cf. Figure 1b).
In all regions, on average, precipitation is varying at different time scales, ranging from seasonal to inter-annual (i.e. 2-8 years; Figure 4b). Continuous wavelet spectra however show that those time scales of variability are non-stationary ( Figure   180 4c), with temporal changes in terms of amplitude discriminating the different regions. For instance, South-western regions are characterized by quasi-continuous annual variability until the late 1980s, while other watersheds show sparsely significant annual variability (Figure 4c). Similarly, although there is significant inter-annual variability in all watersheds from the late and -eastern regions, but are much less pronounced and significant over shorter periods of time than in other regions ( Figure   5).
In summary, different regions with coherent precipitation variability are identified, and are characterised by three time scales of variability: intra-seasonal, annual and inter-annual. The amplitude of those time scales of variability however differs in time and over French territory. The Northern-Alps appear much less coherent, and seems to share spectral characteristics with other 190 regions, especially those characterized by the presence of mountain ranges. Phase-amplitude interactions are presented in Figure 6b. The lower half of the graph, which refers to lower-frequency driving higher-frequency variability, is populated by 5 − 8yr → 2 − 4yr interactions for western and North-centre regions ( Figure   205 6b, pink, yellow, green, red). Centre-eastern regions are also showing lower-frequency variability driving higher-frequency variability, between 8yr → 6yr variability (Figure 6b). Notably, the North-western region is the only one with cross-scale interactions driving the annual cycle (Figure 6a, green).In the upper half of the graph, which refers to higher-frequency driving lower-frequency variability, we only find North-centre and North-eastern regions, showing 2−4 → 4yr and 3−4yr → 7−8yr phase-amplitude interactions (Figure 6b, blue, red). Note that North-centre, North-eastern and Centre-eastern regions show 210 phase-amplitude and phase-phase interactions at very similar time scales (Figure 6a,b, red, blue, black), while time scales of phase-amplitude and phase-phase interactions do not match in Centre-western and North-western and South-western regions ( Figure 6a, b, pink, green, yellow). Regions to the East thus appear to have both phase-phase and phase-amplitude interactions at the same time scales, while western regions are more characterized by phase-amplitude interactions.

cross-scale interactions
The precipitation cross-scale interactions can be of different forms: phase-phase, phase-amplitude, uni-or bi-directional, 215 from lower to higher time scales 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. Inter-estingly, cross-scale interactions tend to converge toward specific time scales, 2-4yr and 5-8yr, which were linked to oceanatmosphere 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). However, the presence of mirror interactions also indicate strong 220 bidirectional negative feedback, which might refer to interaction between ocean-atmospheric variability and land-surface processes (Mariotti et al., 2018). Soil moisture has been shown to be either a positive or negative feedback while several negative feedbacks have been identified in both ocean-atmosphere and land-ocean interactions (Boé, 2013;Sejas et al., 2014).

225
In temperature, nine regions with homogeneous time-frequency patterns are identified ( Figure 7a): North-western-high (pink), North-western-low (black), North-eastern (blue), Centre-eastern (red), Centre-western (green), South-eastern-high (yellow), South-eastern-low (brown), South-western-high (dark green) and South-western-low (purple). Fuzzy clustering shows that watersheds typically converge toward singular clusters, defining highly coherent regions (Figure 7a). This is however not true for the Centre-western region, which is characterized by a mix of spectral characteristic defining other regions (cf. red, green,

cross-scale interactions
240 Figure 9 shows cross-scale interactions for each cluster of temperature variability identified in Figure 8a.
For temperature, phase-phase interactions are mostly concentrated in the upper half of the graph, which refers to higherfrequency modulating lower-frequency ( Figure 9a). Notably, a 2 − 6yr → 6 − 8yr phase-phase interaction appears more pronounced over Northern regions (Figure 9a, blue, red, pink, black). The Centre-western region shows similar phase-phase interactions, but at 1 − 3yr → 4 − 6yr time scales (Figure 9a, green). In the lower half of the graph, which refers to lower-245 frequency modulating higher-frequency, interactions are found at very similar timescales, but at slightly higher frequency, for all regions (e.g. , 2 − 5yr → 1 − 4yr variability). Temperature in the South-western-low region, however, show slightly different characteristics, phase-phase interactions between lower-and higher-frequency occurring between 7 − 8yr → 3 − 4yr and 7 − 8yr → 3 − 4yr variability (Figure 9a, purple).
Temperature phase-amplitude interactions are mostly acting on the 3-4yr time scale for all regions (Figure 9b). In particular, in temperature, more pronounced phase-amplitude interactions are found over the South-western-low region (Figure 9b, purple), which consistent with previous studies on phase-amplitude interactions in European temperature (Palus, 2014;Jajcay et al., 2016). Over South-western regions, temperature, however, shows both 3 − 8yr → 3 − 4yr and 2 − 4yr → 4 − 7yr phaseamplitude interactions (Figure 9b, brown, purple). Furthermore, it should be noted that temperature variability interactions occur between very similar time scales over a number of regions (Figure 9b, pink, red, yellow, purple). According to Paluš

255
(2014), interactions between very similar times, or the same times, can only occur if, at least, two processes are present.
As for precipitation, phase-phase and phase-amplitude cross-scale interactions are region-dependent, and can be uni-or bi-lateral, in temperature. 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 time scales involved for phase-phase and phase-amplitude interactions are very similar in precipitation, they differ largely in tem-260 perature ( Figure 9b). This suggests that the processes driving phase-phase and phase-amplitude cross-scale interactions are different in temperature, and also different between temperature and precipitation. As previously suggested for precipitation, such cross-scale interactions might refer to negative feedbacks between ocean, atmosphere and land-processes, which can be noticeably different in temperature. Temperature variability is mostly derived from changes in heat advection from the North-Atlantic oceans over most of the French territory, except for the Mediterranean regions that are strongly influenced by the 265 Mediterranean Sea (Scherrer, 2006;Sodemann and Zubler, 2010). In addition, temperature's feedback processes (i.e. positive or negative) are extremely region-dependent (Lorenz et al., 2016;Levine et al., 2016;Levine, 2019). Nevertheless, such feedback processes could result from regional changes in soil moisture, as well as from heat exchanges between ocean and land, especially over the coastal regions (Lambert et al., 2011;Sejas et al., 2014).

Time-frequency patterns
Six regions with homogeneous time-frequency patterns are identified in discharge (Figure 10a): North-western (black), Northeastern (blue), North-Centre (red), Centre-western (green), South-eastern (yellow) and South-western (pink). However, several watersheds, especially in the South, show memberships to multiple regions, suggesting lower spatial coherence in discharge than in precipitation and temperature. Lower spatial coherence, however, could mostly be explained by mixing of solid and 275 liquid precipitation in driving discharge variability in the Alps and from the local heterogeneity of precipitation due to convective dynamics in the Pyrenees. 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 time scales, as determined through the global wavelet spectra ( Figure 10b). Unlike other regions, South-eastern watersheds, however, shows significant intra-seasonal variability, even 280 using monthly data (Figure 10b). Continuous wavelet spectra show that both annual and intra-seasonal variability can be non-stationary, with temporal changes in terms of amplitude discriminating the different homogeneous regions of discharge variability (Figure 10c). For instance, annual variability is only significant for specific periods in the South-eastern watersheds, while other regions show mostly constant annual variability (Figure 10c). Similarly, in the South-eastern region, intra-seasonal discharge variability sparsely appears significant from the 1980s, absent from the other regions (Figure 10b).

285
After removing the seasonality, focusing on inter-annual variability, North-eastern watersheds stand out as having continuous significant inter-annual variability throughout the time series, with 4-5yr and 5-8yr variability before and after the 1990s, respectively (Figure 11b). South-eastern and -western regions also stand out, as they show 2-4yr variability in the mid-1970s and 2000s (Figure 11b, yellow, pink). In addition, South-eastern regions do not show significant variability in discharge at time scales greater than 4yr (Figure 11a-b).

290
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 by the presence of mountain ranges (cf. Figure 1b). However, compared to precipitation and temperature, 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.

cross-scale interactions
A major 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 North-eastern, South-eastern and North-western regions (blue, yellow and black; Figure 6a, Figure 12a).
Phase-phase interactions that were identified in temperature are much less evident (Figure 9a, Figure 12a). It should also be 300 noted that many small patches, describing phase-phase interactions in precipitation and temperature, are systematically not transferred to discharge variability (Figures 6a, 9a, 12a). Instead, discharge variability seems to exclusively preserve large patches of phase-phase interactions (Figures 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 the North-centre, where phase-phase interactions are absent in discharge (Figure 12a), but were identified in precipitation and 305 temperature (Figure 6a, 9a).
But, more importantly, there is no phase-amplitude interaction in discharge (Figure 12b). This points out that watershed properties modulate the interacting processes in precipitation and temperature. Since phase-amplitude interactions require a negative feedback system, any change of the feedback sign (by either runoff or subsurface processes) can cancel those phase-amplitude interactions. As our data set is mostly composed of low groundwater support, those modulations are unlikely 310 to results from the water table, especially as phase-phase interactions are inherited from precipitation. In addition, further analysis on Paris' Austerlitz gauging station, which includes very large groundwater support, reveals the same absence of phaseamplitude 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 of phase-320 amplitude interactions.

Discussions an conclusions
As recommended by Blöschl et al. (2019), studying 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 unravelled homogeneous regions of hydroclimate variability in France, accounting for non-stationarity and non-linearity, bringing additional information over 325 previous, regime-based, classifications in France or elsewhere (Champeaux and Tamburini, 1996;Bower and Hannah, 2002;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 spatio-temporal characteristics of each homogeneous region, including characteristic time scales of hydroclimate variability (i.e. precipitation, temperature and discharge) and cross-scale interactions.

330
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 large 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 Tamburini, 1996;Sauquet et al., 335 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 time scales of hydroclimate variability differentiate the regions, highlighting the need for accounting for non-stationary behaviours in global to regional hydroclimate study. Overall, hydroclimate variability displays intra-seasonal (<1yr), annual ( 1yr) and inter-annual (2-4yr and 5-8yr) time scales. Our results, which were focused on the French territory, are therefore consistent with time scales of variability identified over the world major rivers (Labat, 2006).

340
Looking at the internal dynamics of each coherent region, i.e. phase-phase and phase-amplitude causal relationships, complex interactions have been identified. Those interactions can be orientated from lower (higher) to higher (lower) time scales of variability, can be uni-or bi-directional and even self-interacting. For all regions, even regions with very similar time-frequency patterns, the cross-scale interactions are very different, which implies different internal dynamic (Paluš, 2014). As suggested earlier, this might also highlight different interactions between ocean, atmosphere and land-surface processes over the dif-345 ferent regions of France. For instance, the interactions between atmosphere variability (either precipitation or temperature) are known to be strongly modulated by interactions with land-surface processes, such as soil moisture, evapotranspiration, and surface albedo (including snow and ice cover), and these processes are more or less important over the French territory (Ducharne et al., 2020). But more importantly, phase-phase interactions in discharge variability appear partly inherited from precipitation and temperature. This suggests that the watershed characteristics do not entirely filter the interactions between 350 ocean, atmosphere and land-processes driving precipitation and temperature. Interestingly, however, discharge variability does not show any phase-amplitude causality, while such interactions are significant in precipitation and temperature. This suggests that catchment properties substantially modulate the amplitude of the input climate signal (i.e. precipitation and temperature), which are then not transferred to discharge variability. Such finding is consistent with previous studies, showing that infiltration processes and ground water variability acts as a low-pass filter on the input precipitation and temperature variability (Massei 355 et al., 2007). However, our hydrological data set contains very few watersheds with groundwater support. The processes driving those amplitude modulations could not be fully explained within the scope of this study, and require further explorations.
Nevertheless, those findings allow for a better identification of climate deterministic processes controlling hydroclimate variability, notably using cross-scale analysis, which could help 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 360 seamless hydrological predictions based on climate information, as only some parts of the climate signals are transferred to discharge systems. Thus, causal cross-scale relationship could be used to inform and improve existing seasonal to multiyear seamless forecasting for hydrological variability, including extremes (e.g. flood and drought). Preliminary work in this direction were recently proposed by Jajcay et al. (2018), who developed a composite binning method enabling to forecast a particular time series based on conditional phase of another. Similarly, it would be of crucial importance to determine whether 365 hydrological models, which are commonly used in climate-impact assessments, are reproducing the filtering-processes induced by the catchment properties, and identify those (Ducharne et al., 2020). Long term hydroclimate variability only represents a fraction of the total variability, however, strong teleconnections between long term and smaller than inter-annual time scales have been highlighted. Those teleconnections result from both spatial and temporal interactions and feedback between the hydroclimate processes (Feliks et al., 2016). Owing to the recent addition of long term, high spatial resolution hydroclimate 370 data sets (e.g. Fyre reconstructions, Devers et al. (2020Devers et al. ( , 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 time scales. The methodology presented in this work can enable deeper analyses than those based on correlations, which may overlook some important hydroclimate processes.