Hydrol. Earth Syst. Sci., 23, 73–91, 2019
https://doi.org/10.5194/hess-23-73-2019
© Author(s) 2019. This work is distributed under
the Creative Commons Attribution 4.0 License.
A large sample analysis of European rivers on seasonal river flow
correlation and its physical drivers
Theano Iliopoulou
1
, Cristina Aguilar
2
, Berit Arheimer
3
, María Bermúdez
4
, Nejc Bezak
5
, Andrea Ficchì
6,a
,
Demetris Koutsoyiannis
1
, Juraj Parajka
7
, María José Polo
2
, Guillaume Thirel
8
, and Alberto Montanari
9
1
Department of Water Resources and Environmental Engineering, School of Civil Engineering,
National Technical University of Athens, Zographou, 15780, Greece
2
Fluvial Dynamics and Hydrology Research Group, Andalusian Institute of Earth System Research,
University of Córdoba, Córdoba, 14071, Spain
3
Swedish Meteorological and Hydrological Institute, 601 76 Norrköping, Sweden
4
Water and Environmental Engineering Group, Department of Civil Engineering,
University of A Coruña, 15071 A Coruña, Spain
5
Faculty of Civil and Geodetic Engineering, University of Ljubljana, Jamova 2, 1000 Ljubljana, Slovenia
6
Department of Geography and Environmental Science, University of Reading, Reading, RG6 6AB, UK
7
Vienna University of Technology, Institute of Hydraulic Engineering and Water Resources Management,
Karlsplatz 13/222, 1040 Vienna, Austria
8
IRSTEA, Hydrology Research Group (HYCAR), 92761, Antony, France
9
Department DICAM, University of Bologna, Bologna, 40136, Italy
a
formerly at: IRSTEA, Hydrology Research Group (HYCAR), 92761, Antony, France
Correspondence: Theano Iliopoulou (anyily@central.ntua.gr)
Received: 15 March 2018 – Discussion started: 3 April 2018
Revised: 16 November 2018 – Accepted: 6 December 2018 – Published: 7 January 2019
Abstract. The geophysical and hydrological processes gov-
erning river flow formation exhibit persistence at several
timescales, which may manifest itself with the presence of
positive seasonal correlation of streamflow at several differ-
ent time lags. We investigate here how persistence propagates
along subsequent seasons and affects low and high flows. We
define the high-flow season (HFS) and the low-flow season
(LFS) as the 3-month and the 1-month periods which usu-
ally exhibit the higher and lower river flows, respectively. A
dataset of 224 rivers from six European countries spanning
more than 50 years of daily flow data is exploited. We com-
pute the lagged seasonal correlation between selected river
flow signatures, in HFS and LFS, and the average river flow
in the antecedent months. Signatures are peak and average
river flow for HFS and LFS, respectively. We investigate the
links between seasonal streamflow correlation and various
physiographic catchment characteristics and hydro-climatic
properties. We find persistence to be more intense for LFS
signatures than HFS. To exploit the seasonal correlation in
the frequency estimation of high and low flows, we fit a bi-
variate meta-Gaussian probability distribution to the selected
flow signatures and average flow in the antecedent months
in order to condition the distribution of high and low flows
in the HFS and LFS, respectively, upon river flow observa-
tions in the previous months. The benefit of the suggested
methodology is demonstrated by updating the frequency dis-
tribution of high and low flows one season in advance in a
real-world case. Our findings suggest that there is a traceable
physical basis for river memory which, in turn, can be sta-
tistically assimilated into high- and low-flow frequency es-
timation to reduce uncertainty and improve predictions for
technical purposes.
Published by Copernicus Publications on behalf of the European Geosciences Union.
74 T. Iliopoulou et al.: A large sample analysis of European rivers
1 Introduction
Recent analyses for the Po River and the Danube River high-
lighted that catchments may exhibit significant correlation
between peak river flows and average flows in the previ-
ous months (Aguilar et al., 2017). Such correlation is the
result of the behaviours of the physical processes involved
in the rainfall–runoff transformation that may induce mem-
ory in river flows at several different timescales. The pres-
ence of long-term persistence in streamflow has been known
for a long time, since the pioneering works of Hurst (1951),
and has been actively studied ever since (e.g. Koutsoyian-
nis, 2011; Montanari, 2012; O’Connell et al., 2016 and refer-
ences therein). While a number of seasonal flow forecasting
methods have been explored in the literature (e.g. Bierkens
and van Beek, 2009; Dijk et al., 2013), attempts to explic-
itly exploit streamflow persistence in seasonal forecasting
through information from past flows have been, in general,
limited. Koutsoyiannis et al. (2008) proposed a stochastic ap-
proach to incorporate persistence of past flows into a predic-
tion methodology for monthly average streamflow and found
the method to outperform the historical analogue method (see
also Dimitriadis et al., 2016, for theory and applications of
the latter) and artificial neural network methods in the case
of the Nile River. Similarly, Svensson (2016) assumed that
the standardized anomaly of the most recent month will not
change during future months to derive monthly flow fore-
casts for 1–3 months lead time and found the predictive skill
to be superior to the analogue approach for 93 UK catch-
ments. The above-mentioned persistence approach has also
been used operationally in the production of seasonal stream-
flow forecasts in the UK since 2013, within the framework of
the Hydrological Outlook UK (Prudhomme et al. 2017). A
few other studies have included past flow information in pre-
diction schemes along with teleconnections or other climatic
indices (Piechota et al., 2001; Chiew et al., 2003; Wang et al.,
2009). Recently, it was shown that streamflow persistence,
revealed as seasonal correlation, may also be relevant for pre-
diction of extreme events by allowing one to update the flood
frequency distribution based on river flow observations in the
pre-flood season and reduce its bias and variability (Aguilar
et al., 2017). The above previous studies postulated that sea-
sonal streamflow correlation may be due to the persistence
of the catchments storage and/or the weather, but no attempt
was made to identify the physical drivers.
The present study aims to further inspect seasonal persis-
tence in river flows and its determinants, by referring to a
large sample of catchments in six European countries (Aus-
tria, Sweden, Slovenia, France, Spain, and Italy). We focus
on persistence properties of both high and low flows by in-
vestigating the following research questions: (i) what are the
physical conditions, in terms of catchment properties, i.e. ge-
ology and climate, which may induce seasonal persistence in
river flow, and (ii) can floods and droughts be predicted, in
probabilistic terms, by exploiting the information provided
by average flows in the previous months? These questions
are relevant for gaining a better comprehension of catchment
dynamics and planning mitigation strategies for natural haz-
ards. To reach the above goals, we identify a set of descrip-
tors for catchment behaviours and climate and inspect their
impact on correlation magnitude and predictability of river
flows.
A few studies have analysed physical drivers of streamflow
persistence on annual and deseasonalized monthly and daily
time series (Mudelsee, 2007; Hirpa et al., 2010; Gudmunds-
son et al., 2011; Zhang et al., 2012; Szolgayova et al., 2014;
Markonis et al., 2018), but the topic has been less studied on
intra-annual scales relevant to seasonal forecasting of floods
and droughts.
To demonstrate the high practical relevance of the identi-
fied seasonal correlations we present a technical experiment
for one of the studied rivers (Sect. 7) in which the frequency
distribution of both high and low flows is updated one season
in advance by exploiting real-time information on the state
of the catchment.
2 Methodology
The investigation of the persistence properties of river flows
focuses separately on both high and low discharges and is
articulated in the following steps: (a) identification of the
high- and low-flow seasons, (b) correlation assessment be-
tween the peak flow in the high-flow season (average flow
in the low-flow season) and average flows in the previous
months, (c) analysis of the physical drivers for streamflow
persistence and its predictability through a principal compo-
nent analysis (PCA), and (d) real-time updating of the fre-
quency distribution of high and low flows for a selected case
study with significant seasonal correlation by employing a
meta-Gaussian approach. The above steps are described in
detail in the following sections.
2.1 Season identification
Season identification is performed algorithmically to identify
the high-flow season (HFS) and low-flow season (LFS) for
each river time series. For the estimation of HFS, we employ
an automated method recently proposed by Lee et al. (2015),
which identifies the high-flow season as the 3-month period
centred around the month with the maximum number of oc-
currences of peaks over threshold (POT), with the thresh-
old set to the highest 5 % of the daily flows. To evaluate
the selection of HFS, a metric constructed as the percentage
of annual maximum flows (PAMF) captured in the HFS is
used. The PAMFs are classified in the subjective categories
of “poor” (< 40 %), “low” (40 %–60 %), “medium” (60 %–
80 %), and “high” (> 80 %) values, denoting the probability
that the identified HFS is the dominant high-flow season in
the record. If the identified peak month alone contains more
Hydrol. Earth Syst. Sci., 23, 73–91, 2019 www.hydrol-earth-syst-sci.net/23/73/2019/
T. Iliopoulou et al.: A large sample analysis of European rivers 75
than or equal to 80 % of the annual maxima flows, a unimodal
regime is assumed and the identification procedure is termi-
nated. In all other cases, the method allows for the search of
a second peak month and the identification of a minor HFS,
but we do not further elaborate on this analysis here, because
we are only interested in the most extreme seasons for the
purpose of predicting high and low flows.
The method proposed by Lee et al. (2015) has several ad-
vantages that make it suitable for the purpose of this research.
Most importantly, it is capable of handling conditions of bi-
modality, which is usually a major issue for traditional meth-
ods, e.g. directional statistics (Cunderlik et al., 2004). A po-
tential limitation is the assumption of symmetrical extension
of HFS around the peak month, along with the uniform selec-
tion of its length (3-month period). The degree of subjectiv-
ity in the evaluation of the second HFS is another limitation,
which is not relevant here, as we focus on the main HFS.
The LFS is herein identified as the 1-month period with
the lowest amount of mean monthly flow. An alternative ap-
proach of estimating the relative frequencies of annual min-
ima of monthly flow and selecting the month with the highest
frequency as the LFS is also considered.
2.2 Correlation analysis and physical interpretation
through principal component analysis
2.2.1 Correlation analysis
In the case of HFS, a correlation is sought between the max-
imum daily flow occurring in the HFS period and the mean
flow in the previous months, before the onset of HFS. For
LFS, correlation is computed between the mean flow in the
LFS itself and the mean flow in the previous months. We use
the mean flow in the previous month as a robust proxy of
“storage” in the catchment that is expected to reflect the state
of the catchment, i.e. wetter or drier than usual. Since we are
interested in seasonal persistence, we compute the Pearson’s
correlation coefficient for HFS lag up to 9 months and for
LFS lag up to 11 months.
2.2.2 Analysis of physical drivers
Catchment, geological, and climatic descriptors
An extensive investigation is carried out to identify physical
drivers of seasonal streamflow correlation, in terms of catch-
ment, geological, and climatic descriptors.
As catchment descriptors, we consider the basin area (A),
the baseflow index (BI), the mean specific runoff (SR), the
percentage of basin area covered by lakes (percentage of
lakes – PL) and glaciers (percentage of glaciers – PG), and
altitude as candidates for explanatory variables for stream-
flow correlation.
The area A (km
2
) is primarily investigated, as it is repre-
sentative of the scale of the catchment, under the assumption
that in larger basins the impact of the climatological and geo-
physical processes affecting river flow becomes more signif-
icant and may lead to a magnified seasonal correlation.
The BI is considered based on the assumption that high
groundwater storage may be a potential driver of correla-
tion. BI is calculated from the daily flow series of the rivers
following the hydrograph separation procedure detailed in
Gustard et al. (2008). Flow minima are sampled from non-
overlapping 5-day blocks of the daily flow series, and turn-
ing points in the sequence of minima are sought and identi-
fied when the 90 % value of a certain minimum is smaller or
equal to its adjacent values. Subsequently, linear interpola-
tion is used in between the turning points to obtain the base-
flow hydrograph. The BI is obtained as the ratio of the vol-
ume of water beneath the baseflow separation curve versus
the total volume of water from the observed hydrograph, and
an average value is computed over all the observed hydro-
graphs for a given catchment. A low index is indicative of
an impermeable catchment with rapid response, whereas a
high value suggests high storage capacity and a stable flow
regime.
SR (m
3
s
−1
km
−2
) is computed as the mean daily flow
of the river standardized by the size of its basin area. It
may be an important physical driver, as it is an indica-
tor of the catchment’s wetness. PL (%) and PG (%) are
investigated for the Swedish and Austrian catchments, re-
spectively, as lakes and glaciers are expected to increase
catchment storage thus affecting persistence. Lake cover-
age data are based on cartography and are available from
the Swedish Water Archive (https://www.smhi.se/, last ac-
cess: 1 November 2016), while glacier coverage data are
estimated from the CORINE land cover database (https:
//www.eea.europa.eu/publications/COR0-landcover, last ac-
cess: 6 November 2016).
The effect of catchment altitude is also inspected us-
ing relief maps from the Shuttle Radar Topography Mis-
sion (SRTM) data (http://srtm.csi.cgiar.org/, last access:
28 July 2017). The data are available for the whole globe and
are sampled at 3 arcsec resolution (approximately 90 m). To-
pographic information is available for all catchments located
at latitudes lower than 60
◦
north, while a 1 km resolution dig-
ital elevation model is available for Austria.
As geological descriptors we consider the percentage of
catchment area with the presence of flysch (percentage of fly-
sch – PF) and karstic formations (percentage of karst – PK)
for Austrian and Slovenian catchments, respectively, where
this type of information is available. A subset of Austrian
catchments is characterized by the dominant presence of fly-
sch, a sequence of sedimentary rocks characterized by low
permeability, which is known to generate a very fast flow
response. Karstic catchments, characterized by the irregular
presence of sinkholes and caves, are also known for having
rapid response times and complex behaviour; e.g. initiating
fast preferential groundwater flow and intermittent discharge
via karstic springs (Ravbar, 2013; Cervi et al., 2017). Ge-
ological features are also presumed to be linked to persis-
www.hydrol-earth-syst-sci.net/23/73/2019/ Hydrol. Earth Syst. Sci., 23, 73–91, 2019
76 T. Iliopoulou et al.: A large sample analysis of European rivers
tence properties, because geology is the main control for the
baseflow index across the European continent (Kuentz et al.,
2017). PK (%) and PF (%) are estimated from geological
maps of Slovenia and Austria, respectively.
As climatic descriptors, the mean annual precipitation P
(mm year
−1
) and the mean annual temperature T (
◦
C) are
selected. Corresponding gridded data are retrieved from the
WorldClim database (http://www.worldclim.org/, last access:
20 March 2017) at a spatial resolution of 10 arcminutes (ap-
proximately 18.55 km). We note that low mean temperature
regimes are also associated with snow, the presence of which
is also considered in the interpretation of the results. We also
adopt the De Martonne index (IDM; De Martonne, 1926) as
a climatic descriptor, which is given by IDM = P /(T + 10)
and enables classification of a region into one of the fol-
lowing six climate classes, i.e. arid (IDM ≤ 5), semi-arid
(5 < IDM ≤ 10), dry subhumid (10 < IDM ≤ 20), wet subhu-
mid (20 < IDM ≤ 30), humid (30 < IDM ≤ 60), and very hu-
mid (IDM ≥ 60). Additionally, the Köppen–Geiger climatic
classification (Kottek et al., 2006) of the rivers is assessed.
Principal component analysis
To identify which catchment, physiographic, and climatic
characteristics may explain river memory, we attempt to
regress the seasonal streamflow correlation on the physical
descriptors introduced above. We expect the presence of mul-
ticollinearity among the predictor variables, and therefore
PCA (Pearson, 1901; Hotelling, 1933) was applied to con-
struct uncorrelated explanatory variables. In essence, PCA
is an orthonormal linear transformation of p data variables
into a new coordinate system of q ≤ p uncorrelated variables
(principal components – PCs) ordered by decreasing degree
of variance retained when the original p variables are pro-
jected into them (Jolliffe, 2002). Therefore, the first princi-
pal axis contains the greatest degree of variance in the data,
while the second principal axis is the direction which max-
imizes the variance among all directions orthogonal to the
first principal axis, and each succeeding component in turn
has the highest variance possible while satisfying the condi-
tion of orthogonality to the preceding components. Specifi-
cally, let x be a random vector with mean µ and correlation
matrix 6, and the principal component transformation of x
is then obtained as follows:
y = C
T
x
0
, (1)
where y is the transformed vector whose kth column is the
kth principal component (k = 1, 2, . .. , p), C is the p × p
matrix of the coefficients or loadings for each principal com-
ponent, and x
0
is the standardized x vector. Standardization is
applied in order to avoid the impact of the different variable
units on selecting the direction of maximum variance when
forming the PCs. The y values are the scores of each obser-
vation, i.e. the transformed values of each observation of the
original p variables in the kth principal component direction.
PCA has useful descriptive properties of the underlying
structure of the data. These properties can be efficiently vi-
sualized in the biplot (Gabriel, 1971), which is the combined
plot of the scores of the data for the first two principal com-
ponents along with the relative position of the p variables as
vectors in the two-dimensional space. Herein, the distance
biplot type (Gower and Hand, 1995), which approximates
the Euclidean distances between the observations, is used.
Variable vector coordinates are obtained by the coefficients
of each variable for the first two principal components. After
construction of the PCs, a linear regression model is explored
for the case of HFS and LFS lag-1 correlation.
2.3 Technical experiment: real-time updating of the
frequency distribution of high and low flows
In order to evaluate the usefulness of the information pro-
vided by the 1-month-lag seasonal correlation for flow signa-
tures in HFS and LFS, we perform a real-time updating of the
frequency distribution of high and low flows based on the av-
erage river flow in the previous month. A similar analysis for
the high flows was carried out by Aguilar et al. (2017) for the
Po and Danube Rivers. In principle, this is a data assimila-
tion approach, since real-time information, i.e. observations
of the average river flow, is used in order to update a prob-
abilistic model and inform the forecast of the flow signature
of the upcoming season.
In detail, a bi-variate meta-Gaussian probability distribu-
tion (Kelly and Krzysztofowicz, 1997; Montanari and Brath,
2004) is fitted between the observed flow signatures, i.e. peak
flow in the HFS, Q
P
, average flow in the LFS, Q
L
, and the
average flow in the pre-HFS and LFS months, Q
m
. The peak
HFS flow and the average LFS flow are the dependent vari-
ables and are extracted as the peak river discharge observed
in the previously identified HFS and the average river dis-
charge observed in the previously identified LFS, respec-
tively. The average flow in the month preceding the HFS and
the LFS is the explanatory variable in both cases. In the fol-
lowing, random variables are denoted by an underscore and
their outcomes are written in plain form.
The normal quantile transform (NQT; Kelly and Krzyszto-
fowicz, 1997) is used in order to make the marginal probabil-
ity distribution of dependent and explanatory variables Gaus-
sian. This is achieved as follows: (a) the sample quantiles Q
are sorted in increasing order, e.g. Q
m
1
,Q
m
2
. . . Q
m
n
, (b) the
cumulative frequency, e.g. FQ
m
i
, is computed via a Weibull
plotting position, and (c) the standard normal quantile, e.g.
NQ
m
i
, is obtained as the inverse of the standard normal dis-
tribution for each cumulative frequency, e.g. G
−1
(FQ
m
i
).
Therefore, all sample quantiles are discretely mapped into
the Gaussian domain. To get the inverse transformation for
any normal quantile, we connect the points in the above map-
ping with linear segments. The extreme segments are ex-
tended to allow extrapolation outside the range covered by
the observed sample.
Hydrol. Earth Syst. Sci., 23, 73–91, 2019 www.hydrol-earth-syst-sci.net/23/73/2019/
T. Iliopoulou et al.: A large sample analysis of European rivers 77
In the Gaussian domain, a bi-variate Gaussian distribu-
tion is fitted between the random explanatory variable NQ
m
and the dependent variables NQ
P
and NQ
L
by assuming the
stationarity and ergodicity of the variables. We define the
generic random variable NQ
fs
to represent any dependent
flow signature, i.e.; NQ
P
and NQ
L
in our case. Then, the pre-
dicted signature at time t can be written as
NQ
fs
(t) = ρ(NQ
m
,NQ
fs
)NQ
m
(t − h) + Nε(t), (2)
where ρ(NQ
m
, NQ
fs
) is the Pearson’s cross-correlation coef-
ficient between NQ
m
and NQ
fs
, h is the selected correlation
lag with h = 1 in the present application, and Nε(t) is an
outcome of the stochastic process Nε, which is independent,
homoscedastic, stochastically independent of NQ
m
, and nor-
mally distributed with zero mean and variance 1 − ρ
2
(NQ
m
,
NQ
fs
). Then, the joint bi-variate Gaussian probability dis-
tribution function is defined by the mean (µ(NQ
m
) = 0
and µ(NQ
fs
) = 0), the standard deviation (σ (NQ
m
) = 1 and
σ (NQ
fs
) = 1) of the standardized normalized series, and the
Pearson’s cross-correlation coefficient between the normal-
ized series, ρ(NQ
m
, NQ
fs
). From the Gaussian bi-variate
probability properties, it follows that for any observed
NQ
m
(t − h) the probability distribution function of NQ
fs
(t)
conditioned on NQ
m
is Gaussian, with parameters given by
µ(NQ
fs
(t)) = ρ(NQ
m
,NQ
fs
)NQ
m
(t − h), (3)
σ (NQ
fs
(t)) = (1 − ρ
2
(NQ
m
,NQ
fs
))
0.5
. (4)
To derive the probability distribution of Q
fs
(t) conditioned
to the observed Q
m
(t − h), we first apply the inverse NQT,
i.e. we use linear segments to connect the points of the pre-
vious discrete quantile mapping of the original quantiles into
the Gaussian domain, and accordingly, obtain Q
fs
(t) for any
NQ
fs
(t). Subsequently, we estimate the parameters of an as-
signed probability distribution for the obtained quantiles in
the untransformed domain. This is referred to as the up-
dated probability distribution of the considered flow signa-
ture (NQ
P
and NQ
L
, in our case). We use the extreme value
type I distribution for the peak flows and calculate the differ-
ences in the magnitude of estimated maxima for a given re-
turn period between the unconditioned and the updated distri-
bution. The latter is conditioned by the 95 % sample quantile
of the observed mean flow in the previous month. To model
the low flows we use the log-normal distribution, which was
found to exhibit the best fit for the river in question among
other typical candidates for average flows, i.e. the Weibull
and Gamma distribution. The low flows are conditioned by
the lower 5 % sample quantile of the observed mean flow in
the previous month.
3 Data and catchment description
The dataset includes 224 records spanning more than
50 years of daily river flow observations from gauging sta-
tions, mostly from non-regulated streams. A few catchments
are impacted by regulation. Among the 224 rivers, 108 are
located in Austria, 69 in Sweden, 31 in Slovenia, 13 in
France, two in Spain, and one in Italy. Catchment areas vary
significantly, the largest being the Po River basin in Italy
(70 091 km
2
) and the smallest being the Hallabäcken River
basin in Sweden (4.7 km
2
). The geographical location of the
river gauge stations as well as their climatic classification are
shown in Fig. 1. Most of the examined rivers belong to either
a warm temperate (C) or a boreal or snow climate (D) with
a subset impacted by polar climatic conditions (E), accord-
ing to the updated world map of the Köppen–Geiger climate
classification (Fig. 1) based on gridded temperature and pre-
cipitation data for the period 1951–2000 (Kottek et al., 2006).
More specifically, the majority of French and Slovenian and
approximately one third of the Swedish basins belong to the
warm temperate Cfb category characterized by precipitation
distributed throughout the year (fully humid) and warm sum-
mers. The rest of the Swedish catchments are impacted by a
Dfc climatic type, i.e. a snow climate, fully humid with cool
summers. The Austrian catchments belonging to the region
impacted by the European Alps have the most complicated
regime due to their topographic variability. At the lowest al-
titudes, Cfb is the prevailing regime, but as proximity to the
Alps increases, a Dfc regime dominates, and progressively, in
the highest altitude basins, the climate becomes a polar tun-
dra type (Et), characterized primarily by the very low temper-
atures present. The characteristics of all the climatic regimes
of the studied rivers are given in the legend of Fig. 1. A sum-
mary of the river basins under study, in terms of the selected
descriptors, is also provided in Table 1, showing that the in-
vestigated rivers cover a wide range of catchment area sizes,
flow regimes, and climatic conditions.
It is relevant to note that 16 of the Austrian rivers are sub-
ject to regulation, which may alter the persistence proper-
ties of river flows. This relates to generally “mild” forms of
regulation, i.e. upstream regulation with a very low degree
of flow attenuation, hydropower operations, and flow diver-
sions to and from the basin. A preliminary examination of
these rivers did not reveal any significant change during time
of the flow regime. The presence of regulation does not pre-
clude the exploitation of correlation for predicting river flows
in probabilistic terms, but it may affect the analysis of phys-
ical drivers, as it may enhance or reduce persistence in the
natural river flow regime. Given that detailed information is
generally lacking on the impact of regulation (Kuentz et al.
2017), we assume stationarity of the river flows for all the
catchments herein considered and, additionally, assume that
river management does not significantly affect the identifica-
tion of the physical drivers.
www.hydrol-earth-syst-sci.net/23/73/2019/ Hydrol. Earth Syst. Sci., 23, 73–91, 2019
