A novel algorithmic framework for identifying changing streamflow regimes: Application to Canadian natural streams (1966-2010)

Climate change significantly affects natural streamflow regime. To assess alterations in streamflow regime, typically few streamflow characteristics are considered and their significant variations in time and space are taken as a notion of change. Although, this approach is informative, intuitively appealing and widely-implemented, (1) it cannot see simultaneous changes in multiple streamflow characteristics; (2) it does not utilize all the available information contained in a streamflow hydrograph; 10 and (3) it cannot describe how and to what extent one streamflow regime evolves into other regime types. To address these gaps, we conceptualize streamflow regimes as intersecting spectrums that are formed by multiple streamflow characteristics. Accordingly, we recognize that changes in streamflow regime should be diagnosed through gradual, yet continuous changes in an ensemble of streamflow characteristics. To incorporate these key considerations, we propose a fuzzy clustering-based approach to classify the natural streamflow into a finite set of intersecting regime types. Accordingly, by analyzing how the 15 degrees of membership to regime types change, we quantify monotonic shifts between regime types in time and space. Our proposed algorithm eliminates the subjectivity in quantifying shift between flow regimes, and can extract valuable knowledge stored in the shape and variability of annual streamflow hydrographs. We apply this approach to the natural streamflow data, obtained from 106 Canadian gauges, during the period of 1966 to 2010. We show that natural streamflow in Canada can be categorized into six regime types, with clear physical and geographical distinctions. Analyses of trends in membership values 20 during the study period show that alterations in natural streamflow regime are vibrant and can be different within and between major Canadian drainage basins. We show that gradual changes in natural streamflow regimes in Canada can be attributed to simultaneous changes in a large number of streamflow characteristics, some of which have been previously unknown or not well-attended. Our study introduces a generic algorithmic framework for identifying changing streamflow regime at regional and global scales, and provides a fresh look at streamflow alterations in Canada, which can be seen as another line of evidence 25 for the complex and multifaceted impacts of climate change on streamflow regime, particularly in cold regions.


Introduction
Natural streamflow characteristics have been critical consideration for ecosystem and human developments around rivers, globally (Hart and Finelli, 1999;Nazemi and Wheater, 2014;Hassanzadeh et al., 2017). For instance, since https://doi.org/10.5194/hess-2020-334 Preprint. Discussion started: 31 August 2020 c Author(s) 2020. CC BY 4.0 License. early settlements, human learned that timing and duration of low flows control river biodiversity, riparian vegetation and water 30 quality (Rolls et al., 2012;Ireson et al., 2015;Knouft and Ficklin, 2017). During the current "Anthropocene", streamflow characteristics are key factors for infrastructure design as well as land use and land management. While some streamflow characteristics reveal potentials for natural resource development, particularly for agriculture and hydropower production (Hamududu and Killingtveit, 2012;Amir Jabbari and Nazemi, 2019;Nazemi et al., 2020), some others determines consequences of important natural disasters such as flood and drought (Poff and Zimmerman, 2010;Arheimer and Lindström, 35 2015; Rolls et al., 2018;Zandmoghaddam et al., 2019). A set of natural streamflow characteristics determining timing, magnitude, seasonality and inter-annual variability in streamflow time series can collectively define the streamflow regime (Poff et al., 1997). Traditionally, streamflow regimes are considered stationary in time (Milly et al. 2005(Milly et al. , 2008(Milly et al. , 2015. However, the looming effects of climate change along with massive human interventions through land and water management have raised fundamental questions against the feasibility of stationarity assumption for streamflow conditions 40 (Döll and Zhang, 2010;Arnell and Gosling, 2013;Wheater, 2015a, 2015b;Döll et al., 2018). The contemporary literature is full of evidences, revealing major alterations in natural streamflow regime in various regions, induced by heightened climate variability and change (Barnett et al., 2005;Stahl et al., 2010;Rood et al., 2016;Blöschl et al., 2017;Hodgkins et al., 2017;Nazemi et al., 2017;Dierauer et al., 2018). Moreover, projections of future streamflow conditions under climate change conditions show significant alterations in streamflow regime globally in years to come (Zhang et al., 2016;45 Asadieh, and Krakauer, 2017;Eisner et al., 2017;Gizaw et al., 2017;Grantham et al., 2018).
Climate change impacts on natural streamflow regime are more severe in higher latitudes such as in Canada (e.g., Nijssen et al., 2001;Déry and Wood, 2005;Hinzman et al., 2005;Leclerc and Ouarda, 2007;IPCC, 2013;DeBeer et al., 2016;Brahney et al., 2017;MacDonald et al., 2018;Islam et al., 2019;Champagne et al., 2020Dierauer et al., 2020, where the rate of warming is twice of the global average (Bush and Lemmen, 2019). Both observed and projected changes in Canadian 50 streamflow characteristics are subject to significant spatial variabilities (e.g., Burn et al., 2010;Buttle et al., 2016;O'Neil et al., 2017). In northern Canada, for instance, an increase in spring runoff is expected; whereas in southern Pacific, decreasing summer flows are projected (Kang et al., 2016;Curry et al., 2019;Islam et al., 2019). These differences are not only between different regions and/or drainage basins, but can be observed within the same basin and/or between two tributaries with relatively close proximity. For instance, there are significant differences between forms of change in streamflow regime 55 between northern and southern parts of the Pacific (Déry et al., 2009;Monk et al., 2011;Kang et al., 2016;Brahney et al., 2017). Similarly, in northern Canada, glacier-fed rivers show increases in summer runoff (Fleming and Clarke, 2003;Stahl and Moore, 2006); whereas other regime types show tendency toward decreasing summer runoff (Fleming and Clarke 2003;Hinzman et al., 2005;Janowicz, 2008Janowicz, , 2011Foy et al., 2011).
Despite the body of knowledge already gathered around assessing the effects of climate change on streamflow regime, 60 there are still rooms for methodological improvements. Most importantly, among many potential flow characteristics that can constitute and describe streamflow regime, often only few have been taken into account, not only in the Canadian context, but also globally (Whitfield and Cannon, 2000;Regonda et al., 2005;Stewart et al., 2005;Maurer et al., 2007;Blöschl et al., https://doi.org/10.5194/hess-2020-334 Preprint. Discussion started: 31 August 2020 c Author(s) 2020. CC BY 4.0 License. constituted by several streamflow characteristics, and therefore changes in streamflow regimes should be also manifested through changes in an ensemble of streamflow characteristics. Second, we recognize that there are soft rather than hard distinctions between streamflow regimes, and regime shifts occur rather gradual than abrupt. These two considerations lead to development of our proposed algorithm to classify the streamflow regime using a fuzzy clustering approach and monitoring the gradual variations in the degree of belongingness to each cluster in time and space. In brief, we select a set of streamflow 100 characteristics (or features) to collectively characterize the streamflow regimesee Sect. 2.2. We then use the Fuzzy C-Means algorithm (FCM), which is a well-known clustering approach, to classify streams into a set of overlapping regime types during a common baseline periodsee Sect. 2.3. We accordingly quantify changes in degrees of association to each regime type during the entire data period using a moving trend analysis. By monitoring the co-occurrence of divergent trends in membership values, the transitions of regime types to one another can be identifiedsee Sect. 2.4. Finally, we monitor the co-evolution of 105 regime shifts with the alterations in streamflow characteristics through a formal dependency analysissee Sect. 2.5. Figure 1 shows the proposed three-step procedure. Below we describe each step of this algorithm in more details. https://doi.org/10.5194/hess-2020-334 Preprint. Discussion started: 31 August 2020 c Author(s) 2020. CC BY 4.0 License.

Feature selection
Indicators of Hydrologic Alterations (IHAs: Richter et al., 1996) are a set of streamflow characteristics that are commonly applied as features to characterize changes in natural streamflow regime (e.g., Hu et al., 2008;Yang et al., 2012;Wang et al., 115 2018). Different set of IHAs can be considered depending on the application in hand. As the application presented in this paper concerns annual streamflow regime in Canada, see Sect. 3, we consider 15 IHAs, including annual mean flow, monthly mean flows as well as timings of the annual low and high flows. At each stream, we use the mean (first moment) and variance (second moment) of these 15 indicators during a multi-year timeframe to come up with the 30 features that together can capture the shape of expected annual hydrograph and its inter-annual variability during the considered timeframe. Table 1 shows the 120 name and notation of streamflow features used in this study, where =1:15 and =1:15 correspond to mean and variance of the 15 IHAs, respectively.

Fuzzy C-means clustering 125
Clustering is the process of arranging data into a finite set of classes, in a way that members in the same class have similar characteristics. The statistical methodologies used for clustering in hydrology are numerous (see Monk et al., 2011;Olden et al., 2012;Ternynck et al., 2016;Tarasova et al., 2019;Wolfe et al., 2019;Brunner et al., 2020) and have been traditionally limited to non-overlapping (i.e. hard) classes. Recent theoretical developments have relaxed this assumption and considered a set of overlapping (i.e. soft) classes, in particular in the form of fuzzy sets. The association to each fuzzy cluster can quantified 130 using a degree of belongingness, also known as the membership value (see Bezdek, 1981;Sikorska et al., 2015, Piniewski, 2017Knoben et al., 2018). The process of FCM for clustering streamflow regime can be summarized as the following: Assume that flow data of streams during a common timeframe with length of years are available. For each stream, first and second moments of n IHAs (here = 15), i.e. = [ ], = [ ]; ∈ {1, … , }, ∈ {1, … , }, can be extracted during the considered timeframe w. Accordingly, the extracted features can be normalized to avoid scale mismatches in the feature matrix: 135 where ̅ = [ ̅ ] and ̅ = [ ̅ ] are the matrices of Normalized Streamflow Features (NSFs). FCM partitions the streams into fuzzy clusters, such that the sum of distances for all streams ∈ {1, … , } between normalized feature vector and cluster centroids is minimized. This is through an iterative procedure, which aim at finding the cluster centroid by minimizing the 140 following objective function: This objective function is subject to the following two constraints: 145 The number of clusters (i.e. regime types) can be chosen as a priori or empirically using validity indices (e.g., Pal and Bezdek, 1995;Halkidi et al., 2001;Srinivas et al., 2008). Here we implement three validity indices of Xie-Beni index (Xie and 155 Beni, 1991;hereafter ,), separation index (Bensaid et al., 1996;hereafter ,), and partition index (Bensaid et al., 1996;hereafter ) to come up with an optimal number of clusters.

Detection of change in streamflow regimes
Categorizing natural streams into c regime types takes place in a baseline timeframe with the length of l years, in which optimal number of clusters, cluster centroids and the initial membership degrees to each regime type are identified. For each stream, 160 the timeframe can be moved year-by-year and the membership values can be recalculated for the new timeframe using Eq. (3).
This results into c time series of membership values at each stream, showing how the association to each regime type evolves at each streamsee Nazemi et al. (2017) and Jaramillo and Nazemi (2018) for more details on moving window methodology. Figure 2 exemplifies this process in a hypothetical case. In order to quantify the gradual change in membership degrees, the Mann-Kendall trend test with the Sen's Slope estimator (Mann, 1945;Sen, 1968;Kendall, 1975;Hamed, 2008;Gocic and 165 Trajkovic, 2013) can be applied to membership series. As the sum of memberships in each timeframe is one, gradual increase in memberships to one cluster should coincide with gradual decrease in the membership of one or more clusters. This transition can be identified by significant negative dependencies between membership degrees to two clusters at each stream. Assuming the pair of clusters and in stream i, 175 the rate of shift from to can be quantified as: where , ( ) and , ( ) are membership degrees to clusters and in stream during the timeframe ; ∈ {1, … , }; is the number of timeframe; ( , ) and ( , ) are the expected memberships; and ,( , ) is the slope of the best fitted line. 180

Attribution of change in memberships to streamflow features
Changes in membership degrees to each regime type can be attributed to changes in streamflow characteristics. Here, we recognize that the existence of significant dependence between membership values and streamflow features can provides a notion of attribution. Accordingly, we use the Kendall's tau (Genest and Favre, 2007;Nazemi and Elshorbagy, 2012) to detect the co-occurrence between changes in memberships and changes in streamflow features. Figure  The higher 2 is, the larger is the association between changes in the degrees of membership and the considered streamflow characteristics. The coefficient of determination in this case can be calculated as:

Case study and data
Canadian streamflow network is distributed across four major ocean-drained basins, namely Pacific, Atlantic, Arctic and Hudson Bay that together cover around 99.7% of Canada's surface area (Natural Resources Canada, 2007). The Pacific basin, spreading south to north from the US border to Yukon in the west coast, drains the total area of around 1 million km 2 into the 205 Pacific Ocean. The basin is divided from the Arctic basin by the Canadian Rockies. The main sub-basins in the Pacific include Fraser, Yukon, Columbia and the Seaboard. In the east coast, the Atlantic basin drains the total area of 1.6 million km 2 to the Atlantic Ocean. This drainage basin includes important water bodies such as the Great Lakes and is mainly dominated by the streamflow regime in the St. Lawrence River and Seaboard, which are significantly larger sub-basins compared to the Saint https://doi.org/10.5194/hess-2020-334 Preprint. Discussion started: 31 August 2020 c Author(s) 2020. CC BY 4.0 License.
John-St. Croix. Towards north, the Arctic basin drains northern parts of Alberta, British Columbia and Saskatchewan, parts of 210 Yukon as well as the Northwest Territories and the Nunavut. Having the drainage area of over 3.5 million km 2 , the Arctic basin includes some of Canada's largest water bodies such as the Slave, Athabasca and Great Bear lakes as well as Peace, Mackenzie, and Lidar rivers. Mackenzie, Peace-Athabasca and Seaboard are the main sub-basins in the Arctic drainage basin.
With the area of 3.8 million km 2 , the Hudson Bay is the largest drainage basin in Canada, covering five provinces from Alberta in the west to Quebec in the east. The basin includes four major sub-basins, namely Western & Northern Hudson Bay,Nelson,215 Northern Ontario, and Northern Quebec. Nelson, Saskatchewan and Churchill rivers are the major river systems in the Hudson Bay (Pearse et al., 1985;Statistics Canada, 2009). For the sake of convenience, in this paper we consider Northern Ontario and Northern Quebec sub-basins as one single assessment unit.
We use the data from Reference Hydrometric Basin Network (RHBN; Water Survey of Canada, 2017, http://www.wsc.ec.gc.ca/) to diagnose simultaneous changes in natural streamflow regimes across Canadian major basins. 220 RHBN includes 782 stations that measures streamflow at unregulated tributaries over Canada, and therefore are particularly suitable to address climate change impacts on natural streamflow regime (Brimley et al., 1999;Harvey et al., 1999;Whitfield et al., 2012;Zadeh et al., 2020). Considering the available data for the period of 1903 to 2015, we searched for the largest subset of stations with longest continuous daily record during a common period with negligible missing data (i.e., less than a month in a typical year). This results into selecting 106 streamflow gauges during the water years of 1966 to 2010 (1 October 225 1965 to 30 September 2010). Table 2 provides the drainage area, number of stations and abbreviation used for each sub-basin.
For more information about the selected stations, please see Tables S1 to S4 in the Supplement. At each stream, we convert the daily discharge data into runoff depth in millimeter per week by dividing the weekly mean discharge to contributing basin area to obtain a set of comparable streamflow series (Whitfield and Cannon, 2000;Déry et al., 2009). Figure 4 shows the distributions of the considered 106 RHBN stations over major drainage basins and sub-basins. 230 Table 2. Main sub-basins of the four Canadian major drainage basins along with their drainage areas, abbreviations and the number of RHBN stations within their territory used in this study.

Results
We apply the proposed framework to 106 selected RHBN streams using the thirty streamflow features introduced in Table 1.
Note that in this section we only consider a decadal timeframe to calculate the streamflow features, cluster centers, and membership values. We repeat the algorithm with 15-and 20-year timeframes and compare the findings with those presented 240 here in Sect. 5.2 to address the sensitivity of our results to the length of timeframe.

Identifying natural streamflow regimes in Canada
By having a decadal timeframe, we consider the period of 1966-1975 as the baseline. The optimal number of clusters is identified empirically from the pool of = {2, 3, …, 10}, using the three validity indices introduced in Section 2.3. Figure 5 summarizes the results, showing = 6 as the optimal number of clusters. Table 3 introduces these six regimes along with their 245 notation and architype streams. Architype streams are those streams that have the highest association to the identified regime types and can represent the characteristics of a given regime better than other members of the cluster.   Figure 6 shows the distribution of streams belong to each flow regimes across Canada, where red stars represent the architype stations for each regime type. The larger the size of circles are, the higher the degrees of membership are to each 255 cluster. Note that for each regime type, we only consider/show streams with membership values of 0.1 and higher to eliminate streams with insignificant memberships and be able to provide a synoptic look at distribution and physical characteristics of streams within each cluster. As Fig. 6 shows the identified clusters are geographically concentrated and referring to the alreadyknown regimes across the country (see Whitfield, 2001;Bawden et al., 2015;Najafi et al., 2017;Bush and Lemmen, 2019). 260 The first cluster (C1) resembles glacial flow regimes, characterized by gradual rise after spring snowmelt, prolonged peak discharge throughout summer, gradual recession during fall, and low runoff in winter (Déry et al., 2009). Glacial regime spreads mostly over Pacific (34%), Arctic (32%) and Hudson Bay (32%). The Kazan River releasing into the Baker Lake in Nunavut is the architype stream for this regime type. The second cluster (C2) resembles nivo-glacial streamflow regime, dominated by combined melts of accumulated snow 270 and glacier storage (Eaton and Moore, 2010). The hydrograph shape is similar to glacial regime, but with an earlier and sharper recession. Peak flow usually occurs in June and July and the glacier melt maintain flows during late summer and early fall (Moore et al., 2012;Schnorbus et al., 2014). Nearly 46% of nivo-glacial streams are located in Pacific, especially in Fraser and Columbia. The rest are distributed in Hudson Bay (nearly 22%), particularly around the Rocky Mountains, Arctic (19%) and Atlantic (13%). The Clearwater River near Clearwater in southern Alberta is the representative stream of nivo-glacial 275 regime among the considered RHBN stations.
The third cluster (C3) features nival regime, identified by short high flow period in early spring, sharp recession in summer and high flow variability throughout a typical year. 90% of nival streams are located in Atlantic, particularly around southwestern part of St. Lawrence. The rest of streams with nival regime located in southern part of Pacific and northern Ontario. The Matawin River originated from the lake Matawin in Quebec is the architype for nival regime. 280 The fourth cluster (C4) represents nivo-pluvial regime, with two distinct peaks, one in fall due to high precipitation, and one in spring induced by snowmelt. In this regime, snowmelt is the main contributor to runoff; and therefore, the maximum annual discharge is in the spring. This regime has a low flow in early fall and late winter (Hock et al., 2005). 80% of streams with nivo-pluvial regime are located in the southern Atlantic, with few streams in the southern Pacific, and in Hudson Bay.
Gander River at Partridgeberry Hill in eastern Newfoundland is the architype for this regime. 285 The fifth cluster (C5) resembles the hybrid pluvio-nival regime, in which annual stream is influenced more by rainfall around later fall, followed by light increase in discharge due to snowmelt in early spring (Kang et al., 2016). The majority of streams with pluvio-nival regime (81%) are located in lower parts of Atlantic Seaboard. The other 19% of streams are located in Pacific, mainly in southern Fraser and Pacific Seaboard. Beaver Bank River in Nova Scotia is the representative stream for this regime type. 290 Finally, the sixth cluster (C6) features the pluvial regime, in which runoff is dominated by heavy precipitation, especially during winter and lower runoff during summer (Wade et al., 2001;Whitfield, 2001). Due to variations in rainfall, there is a high variability in winter flows. Pluvial regime is only seen in the Vancouver Island with the Sproat River near Alberni being the architype stream.

Detection of change in streamflow regimes 295
To highlight the changes occurred in each regime during the study period, we first look at changes in the shapes of annual streamflow hydrographs in the architype streams, revealed in two distinct decadal episodes during the first and the last decadal timeframes (i.e., 1966 to 1975 vs. 2001 to 2010).  To have a better look at shifts in streamflow regime throughout the country, we calculate the decadal membership values 310 in all streams and throughout the study period and apply the Mann-Kendall trend test with the Sen's Slope on the time series of decadal memberships. Figure 8 summarizes the results. In this figure, rows and columns are related to sub-basins (see Table   2) and regime types (see Table 3 are different between the northern (increasing trends) and southern parts (decreasing trends). Such sub-regional differences in the form of change are quite common. For instance, while in Columbia the degrees of membership to nivo-glacial regime mainly decrease, mixed patterns of change are observed in the Fraser.
The same analysis can be applied to other basins. After the Pacific, the Atlantic basin shows the highest diversity in the streamflow regime by having five regime types within its territory. The dominant streamflow regimes in the Atlantic basin are 330 nival and nivo-pluvial regimes, where 46 and 41 streams out of the total of 50 natural streams considered in this basin show some levels of belongingness to these two regime types, respectively. From these 46 and 41 streams, 29 and 19 streams show increasing trends in belongingness to nival and nivo-pluvial regimes, respectively. The mixed pattern of change is subject to large geographic variations at the sub-basin scales. For instance in the Saint John-St. Croix, the associations to nivo-pluvial regime become stronger during the study period. Similarly, half of the 50 considered streams in Atlantic are associated to some 335 extent with pluvio-nival regime. These streams also show a mixed pattern of change in their belongingness to this regime, similar to the glacial and nivo-glacial regimes that are less common in the Atlantic basin.
In contrast to the Pacific and the Atlantic, the Arctic has the least diversity in the streamflow regime. All considered 12 streams are associated with large degrees to glacial regime, out of which five and six streams show increasing and decreasing trends in the membership, respectively. In addition, nine out of 12 streams located in the Arctic are associated with nivo-glacial 340 regime, for which the degrees of membership increase in six stations. There is only one stream in the Peace Athabasca, showing an increase in belongingness to nivo-pluvial regime during the study period.  In the Hudson Bay basin, 12 and 11 streams out of 14 streams show some levels of association with glacial and nivo-glacial regimes, out of which nine and seven streams show departure from these regimes during the study period, respectively. Six out of 14 considered streams in this basin are associated with some degrees to nival and nivo-pluvial regimes. These streams 350 are located mainly in the Western & Northern Hudson Bay and Nelson sub-basins. The membership values for both regime types show increasing trends, in which three streams exhibit significant trends.

Identifying forms of transformation in streamflow regimes
As noted in Sect. 2.4, representing streamflow regimes using fuzzy clusters implies that the sum of memberships to all regimes at each stream and during each timeframe equals to onesee Eq. 2(b). As a result, the existence of trend in membership values 355 can be taken as a direct evidence for transformation of flow regime from one type to others. This sets the scene for investigating regime shifts quantitatively, as the sign and the slope of change in membership degrees can provide a notion for the direction and intensity of the shifts between regime types. Figure 9 summarizes the results of this analysis for the time series of decadal memberships shown in Fig. 8 in grey. Rows display the considered 106 RHBN stations, categorized by their corresponding basins and sub-basins (see Table 2). The streams in each sub-basin are sorted north to south from top to bottom, similar to Fig.  360 8. Columns are related to the six regime types (see Table 3). Similar categorization is made on the x-axes in all panels.
Considering each panel, the regime type identified by the column filled with diagonal lines (i.e., identical regime types) is the potential receiver regime that the increase in its membership should coincide with the decrease in the membership values of dissimilar regime types (i.e., givers), shown in the x-axes. Considering each pairs of dissimilar regime types, shades of grey represents dissimilar directions of change, quantified according to the sign of Kendall's tau. Note that we only display those 365 pairs of flow regimes, in which evolutions in membership degrees are significantly dependent according to Kendall's tau (pvalue ≤ 0.05) for streams with membership degrees of 0.1 and higher. In each panel, cells related to identical regime types are shaded in diagonal lines.
In the Pacific, the dominant regime shift is observed between glacial and nival regime (in 13 streams; panels of P1 to P4 under glacial and nival). At the sub-basin scale, however, the dominant regime shift may be different. In Yukon, gradual shifts 370 are observed from nival to glacial regime (panel P1/glacial); whereas, in the northern Pacific Seaboard, one stream depart from glacial and shift toward nivo-glacial (panel P2/nivo-glacial). The other stream in this region mainly shift from pluvio-nival to nivo-glacial regimes. In the southern part of Pacific Seaboard, a mixed pattern of transition happens mainly between pluvionival and pluvial regime. In Fraser, a mixed pattern of transition between glacial and nivo-glacial happens in the northern part (panel P3/glacial and nivo-glacial); whereas, one stream in the southern part shift from nivo-glacial to pluvio-nival regime 375 (panel P3/pluvio-nival). In Columbia, two streams located in the northern part shift from nivo-glacial to glacial regimes (panel P4/glacial); while, streams located in the southern part show a strong tendency to depart from glacial to nival as well as from glacial to nivo-pluvial (panels P4/nival and nivo-pluvial). Additionally, a coherent pattern of shift from nivo-glacial to nival as well as from nivo-glacial to nivo-pluvial is observed for streams located in Columbia (panels P4/nival and nivo-pluvial).
This has led to disappearance of glacial and nivo-glacial regimes in some of the streams in this regionsee Fig. 8. In the Atlantic, dominant transition is between nivo-pluvial and pluvio-nival (in 27 streams; panels of At1 to At3 under nivo-pluvial and pluvio-nival), as well as between nival and pluvio-nival regimes (in 25 streams; panels of At1 to At3 under nival and pluvio-nival). At the sub-basin scale, however, the pattern of transition between flow regimes is not homogenous. In the northern part of Atlantic Seaboard, a coherent pattern of shift from pluvio-nival to nival regime is observed for five out of 390 10 streams (panel At1/nival); whereas, four streams in this region show an intense shift mainly from pluvio-nival to nivopluvial regime (panel At1/nivo-pluvial). Similar to the streams in the northern part of Atlantic Seaboard, nine streams located in the southern part shift from pluvio-nival to nivo-pluvial. Moreover, 12 streams in this regions shift from pluvio-nival to nival regime. There are also three stations gradually shift from nival to pluvio-nival regime in this region. In the northern part of St. Lawrence, a shift mainly from nival to glacial is seen in two streams. Reversely, in the southern part seven streams shift 395 mainly from glacial to nival regime (panel At2/nival). Additionally, six streams show tendency to shift mainly from glacial to nivo-pluvial regime. In Saint John-St. Croix, similar to the streams located in Atlantic Seaboard, five streams shift from pluvio-nival to nivo-pluvial regime and the rates of transition in northern streams are more pronounced than the southern parts (panel At3/nivo-pluvial).
In the Arctic, there are frequent transformations between glacial and nivo-glacial regimes (eight out of 12 streams; panels 400 of Ar1 to Ar3 under glacial and nivo-glacial regimes). Similar pattern of shift is observed at the sub-basin scale. In the Seaboard and Lower Mackenzie, the dominant form of transition is from glacial to nivo-glacial regime (see panels of Ar1 and Ar2 under glacial and nivo-glacial regimes). In the Peace Athabasca, two streams shift from nivo-glacial to glacial, while the other slightly departs from glacial into nivo-pluvial regime (panel Ar3/nivo-pluvial).
In the Hudson Bay, the dominant transitions are between glacial and nival regimes (seven out of 14 streams; panels H1 to 405 H3 under glacial and nival regimes), as well as between nivo-glacial and nival regimes (seven out of 14 streams; panels H1 to H3 under nivo-glacial and nival regimes). This is consistent with the pattern of shift at the sub-basin scale. In the Western & Northern Hudson Bay, there are obvious shifts mainly from glacial to nival as well as from nivo-glacial to nival (panel H1/nival). Additionally in this region, a shift from glacial to nivo-pluvial as well as from nivo-glacial to nivo-pluvial is observed. In Northern Quebec & Ontario, natural streamflow regime shifts from glacial toward nivo-glacial regime (panel 410 H2/nivo-glaical). In the northern part of Nelson, three streams shift mainly from glacial to nival regime (panel H3/nival).
Unlike the northern part, one stream in the southern part shifts from nivo-glacial to glacial regime. The other stream also show a tendency to shift from nival to glacial regime.

Attribution of regime shift to changes in streamflow characteristics
The methodology presented in Sect. 2.5 provides an empirical basis to evaluate the association of the regime shifts to changes 415 in streamflow characteristics using the coefficient of determination. Figures 10 and 11 below show this analysis between decadal memberships and decadal means and variances of the 15 IHAs, respectively. In both figures, each panel corresponds with a sub-basin and a regime type, arranged in rows and columns, respectively. In each panel, streamflow features are displayed in the x-axis, i.e., the means of the 15 IHAs are displayed in Fig. 10 and variances are shown in Fig. 11. y-axes show the streams at each basin, sorted from north (top)

475
In the Arctic Seaboard, earlier but more variable timing of annual high flow, increasing mean and variability of seasonal flow in fall, along with increasing winter flows and heightened variability in monthly flow in June collectively correspond with the regime shift from glacial into nivo-glacial regime (see panels Ar1/glacial in Figs. 10 and 11). In Lower Mackenzie, increasing annual and seasonal flows during fall and winter, increasing June's monthly flow, along with heightening variability in the timing of high flows correspond with a gradual shift from glacial to nivo-glacial regime (see panels under Ar2/glacial 480 in Figs. 10 and 11). Having said that, some other streams in this region have experienced a reverse shift from nivo-glacial to https://doi.org/10.5194/hess-2020-334 Preprint. Discussion started: 31 August 2020 c Author(s) 2020. CC BY 4.0 License. glacial regime, due to decreasing monthly flows in October and August as well as decreasing annual flow. In Peace Athabasca, the reverse shift from nivo-glacial to glacial regime corresponds with the earlier and less variable timing of low flows as well as decreasing monthly flow in July (see panels under Ar3/nivo-glacial in Figs. 10 and 11). The other stream in this region shifts from glacial to nivo-pluvial regimes, due to increasing monthly flow in June, decreasing monthly flow in May, increasing 485 mean and variation of monthly flow in March along with delayed and more variation in timing of low flow.

Summary of findings and positioning against earlier studies
By having four ocean-drained basins, spread over more than 6% of the global land area, Canada exhibits a large diversity in its natural streamflow regime. Canadian rivers, rolling coast to coast to coast, support important socio-economic activities such 500 as agricultural and hydropower production, and feed some of the world's most important freshwater bodies that are home to various (some endangered) wildlife species. Natural streamflow regime in Canada, however, is going through drastic changes in recent years. The literature of Canadian hydrology is rich in terms of documenting changes in streamflow characteristics across the country. There are a large body of studies, reporting shift in streamflow regime across different regions due to changes in temperature patterns, magnitude and form of precipitation, snowmelt and snow accumulation processes as well as 505 glacier retreat and permafrost degradation. Thanks to pioneering works of so many before us, including the iconic late northern hydrologist, Richard Janowicz, to whom this paper is dedicated. Having said that, to the best of our knowledge, our work is the first study in which a fully algorithmic framework is used to provide a temporally-homogeneous pan-Canadian view on the recent shifts in natural streamflow regime across the country.
Our results presented in Sect. 4.4 reveal that shifts in streamflow regimes can be attributed to simultaneous changes in a 510 large ensemble of streamflow characteristics. This conclusion is consistent with the earlier findings on changes in natural streamflow characteristics across most Canadian basins and sub-basins; yet our study reveals some new changes in streamflow https://doi.org/10.5194/hess-2020-334 Preprint. Discussion started: 31 August 2020 c Author(s) 2020. CC BY 4.0 License. characteristics that have been previously overlooked or remained unknown. To better position our results against earlier studies, Table 4 summarizes our findings in terms of dominant regime shifts and associated changes in streamflow characteristics at the sub-basin scale. Note that in the majority of sub-basins, there is an obvious divide between dominant 515 regime shifts in northern and southern regions; and therefore, they are separated from one another. Table 4 also makes a clear distinctions between the earlier findings that reconfirmed in the current study, and those exclusively found in our work.
Two important findings can be obtained from this comparison. First and foremost, our study provides new sets of understanding on shift in streamflow regime and forms of alteration in streamflow characteristics in two regions in Canada that have not been previously studied, i.e. northern Fraser and southern St. Lawrence. In both regions, shifts in streamflow 520 regime is attributed to changes in multiple streamflow characteristics. Second, earlier studies often looked only at the changes in expected monthly, seasonal and annual flows. In fact, evaluating changes in variability in streamflow characteristics remained mainly limited to timing of low and high flows. Our study clearly shows that changes in variability of monthly, seasonal and annual flows can be important drivers of shift in streamflow regime across majority of sub-basins in Canada. This is another line of evidence for the complex and multifaceted nature of change in streamflow regime, and the need for 525 simultaneous look at alterations in both expected values and variability of streamflow characteristics to diagnose changes in natural streamflow regime. It should be also noted that earlier studies may have different study periods, and may include streams that are not particularly within the RHBN streams.

Addressing uncertainty
The results provided in Sect. 4 are based on decadal timeframes. To address the uncertainty in our findings due to this 530 assumption, we repeat implementing the proposed algorithm with 15-and 20-year timeframes, and look at the differences in our findings presented in Sect. 4.1 to 4.4.
In terms of clustering, our results show that the number of optimal clusters does not change by altering the timeframe's length. In addition, as shown in Fig. S1 in the Supplement, there are not significant changes in the cluster centers. This highlights the robustness of the FCM in identifying distinct flow regimes. 535 We also look at possible differences in the direction of trends in membership degrees, presented in Sect. 4.2 (Fig. 8).    (Déry et al., 2009;Pike et al., 2010) Delayed and more variable timing of annual low flow; increasing variability in February's monthly flow Fraser (north) Case 1: glacial to nivo-glacial; Case 2: nivoglacial to glacial No earlier study in this region was found. Increasing annual and winter flows (Smith et al., 2007;Walvoord and Striegl, 2007;St. Jacques and Sauchyn, 2009;Rood et al., 2016) Increasing annual and seasonal flows during fall; increasing June's monthly flow; heightening variability in the timing of high flows Peace Athabasca Nivo-glacial to glacial Decreasing monthly flow in July (Yang et al., 2015) earlier and less variable timing of low flows

Hudson Bay
Western & Northern Hudson Bay Glacial to nival Increasing winter flows; decreasing summer flows; increasing variability in winter flows (Déry et al., 2011(Déry et al., , 2018 Delayed and more variable timing of low flows; increasing variability in February's monthly flow The results are presented in Figs. S4 and S5 in the Supplement and are intercompared with corresponding results obtained by decadal timeframes in Fig. 12 (middle column). Our analysis shows that results obtained by 15-and 20-year timeframes are in 550 large agreements with the results obtained using decadal timeframes. Even for the case with the largest disagreement (i.e. nivopluvial regime in the Atlantic), there is 86% agreement in terms of direction of shift in streamflow regimes, obtained by 10and 20-year timeframes.
In terms of attribution of regime shifts to changes in streamflow characteristics, again the results obtained by 15-and 20year timeframes (see Figs. S6 and S7 as well as S8 and S9 in the Supplement, respectively) show large agreement with the 555 results obtained by the decadal timeframe and presented in Sect. 4.4. According to the intercomparison made in Fig. 12 (right column), there is at least 80% agreement in the results obtained by different timeframe lengths. https://doi.org/10.5194/hess-2020-334 Preprint. Discussion started: 31 August 2020 c Author(s) 2020. CC BY 4.0 License.

Conclusions and outlook
This study presents an attempts toward a generic algorithmic framework for identifying changing streamflow regimes, not only in Canada but also globally. The proposed approach is based on two fundamental considerations. First, we recognize that 565 streamflow regime is collectively formed by a large set of streamflow characteristics that can describe the expected annual streamflow hydrograph and the inter-annual variability around it. Second, we acknowledge that streamflow types are rather spectrums, not clear-cut states; and the transition from one regime type to another is gradual rather than abrupt. To accommodate these two considerations, we suggest representing streamflow regime types as intersecting fuzzy sets, in a way that the belongingness of each stream to each regime type can be quantified by a unique membership function. Accordingly, 570 monitoring the trends in membership values in time and space can provide a basis to identify the regime shift from one type to another. In addition, analyzing the covariance of membership values with streamflow characteristics can provide a basis to attribute the regime shift to alterations in certain streamflow characteristics in time and/or space.
By applying the proposed procedure to 106 RHBN streamflow gauges in Canada, we provide a comprehensive look at forms and extents of change in natural streamflow regime during 1966 to 2010 throughout the country. We show that the 575 considered natural streams in Canada can be categorized into six distinct regime types with clear physical and geographical interpretations. Analyses of trends in membership values show that alterations in natural streamflow regime are vibrant and can be different across major drainage basins in Canada. Overall, dominant modes of transition at the basin scale are between glacial and nival in the Pacific, between nivo-pluvial and pluvio-nival as well as nival and pluvio-nival in the Atlantic, between glacial and nivo-glacial in the Arctic, and between glacial and nival as well as nivo-glacial and nival regimes in the Hudson 580 Bay. The details of change in streamflow regime, however, are subject to a large spatial variability within each drainage basin.
For instance in the Pacific, the association to the glacial regime is increasing in Yukon and northern parts of Columbia and Fraser sub-basins; but it is significantly decreasing in the southern regions. This can be due to different manifestations of climate change, which are more revealed as temperature increases in the northand therefore more glacial retreatrather than growing ratios of rain over precipitation in the south, which shift the streamflow more toward rain-dominated regimes 585 (Fleming and Clarke, 2003). Such north/south divides are observed in other drainage basins as well, e.g. in the Atlantic and between streams located in north and south of Bay of Fundy. Having said that, even within close proximity, e.g. between Yukon and northern Pacific Seaboard, differences in the evolution of streamflow regime can be significant. It has been pointed that this regional contrast can be due to existence of the larger glacier in Yukon, i.e., St. Elias Mountains, which exhibit a different response under the same climatic conditions (Fleming and Clarke, 2003;Stahl and Moore, 2006). This reconfirms 590 the important role of landscape in regulating the streamflow response to climate change.
Climate-driven changes in the streamflow regime will have multiple impacts on socio-economic activities and ecosystem services in Canada and globally, and therefore should be managed with care. This requires understanding that water does not respect political boundaries. We strongly believe that Canada and the rest of the world need sustainable and well-integrated water management plans to face the challenges (and the opportunities) ahead of us during the current "Anthropocene". 595