Introduction

Rivers are heterogeneous landscapes (Wiens 2002; Ward et al. 2002; Gilvear et al. 2016). Much like in terrestrial landscapes, the heterogeneity of rivers implies that they are characterised by distinct patches separated by boundaries (Naiman et al. 1988; Poole. 2002; Thorp et al. 2008), as opposed to continuous longitudinal systems (cf. Vannote et al., 1980). Longitudinal boundaries within river networks occur when the structural and or functional properties of adjacent river reaches change discontinuously or non-monotonically in space and time (Yarrow and Marin, 2007). Boundaries are generally not immediate transitions between reaches (or ‘patches’) of a river network but ‘critical zones’ of transition from the conditions within one reach to those of another (Delcourt and Delcourt 1992; Forman 1995). Boundaries influence flows of energy, materials, and organisms between reaches and through a river network (cf., Hitt and Angermeier 2008). They can elicit abrupt or gradual changes in river network character between adjacent reaches (Thorp et al. 2008) and often are associated with increased morphological heterogeneity and species diversity compared to the adjacent habitat patches themselves (Forman 1995). It has also been suggested that boundaries are important for the stability of networks (Stewart 2004), resilience to disturbance (Ash and Newth 2007; Rodriguez-Iturbe et al. 2009), maintaining biodiversity (Grant et al. 2007), and as indicators of environmental change (Naiman et al. 1988) within river networks. Boundaries are therefore a key feature influencing the riverine landscape continuum and, conversely, discontinuum.

Boundaries in river networks can be anthropogenic or natural. The former includes structures constructed across river systems (e.g., dams, weirs and road crossings) with their locations being site specific in terms of human activity. Natural boundaries on the other hand exist in a wide range of environments within river networks and are commonly associated with tributary junctions (river confluences), lithological controls such as waterfalls, and where local geomorphological conditions contribute to significant and often abrupt changes in downstream hydrological and sediment processes, and morphological character. The presence of boundaries within a river network contributes to a longitudinal discontinuum—anthropogenically, with the effects of dams (Stanford and Ward 2001)—and naturally, through river zones with different geomorphic structures and process drivers (Poole 2002). In terms of the natural discontinuum, a dominant paradigm within river science is that the presence of tributaries, albeit of different relative size, determines the character of the river discontinuum. The Network Dynamics Hypothesis (NDH) of Benda et al. (2004) eloquently describes how tributary confluence effects may vary in terms of the specific attributes of a network’s structure. The basic thesis of the NDH is that the probability of significant morphological change in main stem channels increases with the ratio of tributary size to main stem channel size. Accompanying the development of the NDH is a series of testable predictions, most of which have not been thoroughly assessed.

Confluences are not the only type of natural boundaries to occur within river networks. Riverine landscapes are increasingly being viewed as compositions of hierarchically nested patches displaying a high degree of internal heterogeneity in space and time (Petts and Amoros 1996; Montgomery 1999; Thorp et al. 2008). At the drainage basin scale (> 102 km), riverine landscapes have been quantitatively characterized as series of distinct river zones—functional process zones (FPZs) that are large tracts of the river network with similar hydro-geomorphological character throughout a river network (Thoms et al. 2018). Critical transition zones exist between adjacent FPZs as a result of local and regional influences. Phillips (2008) describes five key transition zones between the six main river zones of the Lower Sabine River, USA. These transition zones represent significant boundaries between the river zones in terms of their geomorphology, the majority of which did not occur at tributary junctions. Controls on these zones were inferred to be the result of static (lithological), dynamic (changes over time), and or chronic or continuous influences (Phillips 2008). River networks do exhibit emergent properties, whereby structure and function at higher levels (e.g., entire watersheds) cannot be simply deduced from the collective knowledge of their parts at lower levels (e.g., individual reaches and or tributary junctions) (Allen and Starr 1982).

Quantitative studies of boundaries in riverine landscapes have, until recently, focused on tributary junctions (cf. Benda et al. 2004; Rice et al. 2008) or artificial boundaries such as dams (Stanford and Ward 2001). Less attention has been paid to geomorphic boundaries at other locations. Boundaries or discontinuities between river zones in river networks have been recognised as being important conceptually (Poole 2002) and the multi-scale implications of boundaries to river ecosystems have been addressed by Bretschko (1995) and Ward et al. (1998); however, few have quantified their distribution or strength of influence on the river discontinuum. Boundaries occurring at tributary confluences have been investigated under a stream ordering paradigm of downstream change (Rhoads 1987; Benda et al. 2004); however, many factors influence their occurrence throughout river networks, not just tributary confluences and stream order (cf. Poole, 2002; Thorp et al. 2008). The focus so far on tributary junctions and artificial boundaries in river networks obscures other potentially important boundaries that are likely important for understanding and managing connectivity and fragmentation in river networks.

This study investigates the spatial organisation and strength of boundaries between geomorphologically distinct patches (FPZs) within the river networks of 10 drainage basins in the Kimberley region of north-western Australia. Functional process zones were quantitatively characterised and the boundaries between these geomorphologically distinct zones identified at the river network scale. Boundary strength was determined statistically and based on the geomorphic contrast between the adjacent FPZs, enabling the identification of possible areas of increased physical habitat diversity, as hypothesised by Benda et al. (2004). Possible drivers of the spatial organisation of boundaries throughout the river networks of the Kimberley region, including tributary junctions and lithologies, are also analysed and discussed.

Study area

The Kimberley region, located in north-western Australia (Fig. 1), is comprised of ten drainage basins (or drainage regions, in the case of Cape Leveque), which flow to the Timor Sea and Indian Ocean. These ten basins range in area from 9,631 to 95,344 km2 with a combined catchment area of approximately 306,100 km2. The river network of these basins has a total length of 35,746 km, at the 1:250,000 scale, with drainage densities ranging between 0.03 and 0.14 km per km2 for the individual basins. These basins are amongst the least disturbed by European occupation in Australia (Stein et al. 2002) and are considered to be in a relatively pristine state (Halse et al. 2002).

Fig. 1
figure 1

The Kimberley region of Australia showing the ten major drainage basins in the region, selected place names, and mountain ranges

Full size image

Five broad lithologies exist across the Kimberley region and these are sedimentary rocks; granites; mafic and felsic volcanics; granulite-facies metamorphics; and a mix of mafic-ultramafic intrusives, dolerites and gabbros. Sedimentary rocks cover > 76% of the region, with the remaining four lithological groups contributing varying amounts to the individual basins. Some of the lithological formations in the Kimberly region have been dated at 3.5 billion years before present and this antiquity is because of the region’s tectonic stability (Petheram and Kok 1983). Examples include the Archaean metamorphosed sandstones, quartzites and schists of the Carr Boyd Range in the eastern Kimberley, as well as some exposed Archaean granites on the south-eastern fringes of the King Leopold Range (Petheram and Kok 1983).

The Kimberley region has a tropical monsoonal climate with distinct ‘wet’ (November – March) and ‘dry’ (April – October) seasons. Approximately 90% of annual rainfall (regional long-term mean annual rainfall = 979.2 mm) occurs during the ‘wet’ season when the Kimberly region is influenced by tropical cyclones and low pressures systems. A distinct NW – SE rainfall gradient exists across the region; varying from a long-term annual mean of 1500 mm near Kalumburu, to 400 mm near Broome. Mean monthly maximum temperatures are spatially uniform throughout much of the Kimberley, ranging from 30 oC in the ‘dry’ season to 38 oC in the ‘wet’ season. Local topographic variations and the proximity to the coast can contribute to some variation.

Flow regimes across the Kimberly region reflect the highly variable rainfall patterns. Typically, flows are highly intermittent and most rivers experience between 100 and 200 days per year of no flow, on average; although some rivers experience > 250 no flow days per year. Marked spatial differences in river flow regimes also occur. The Fitzroy River, for example, has a mean monthly flow of < 1,000 ML/day in its headwater reaches compared to month flows > 3,800,000 ML/day in its lower reaches during the ‘wet’ season (Department of Water 2010). There are two dominant flow regime classes for rivers in the Kimberley Region (Pusey et al. 2009). The predictable summer highly intermittent flow regime rivers that exhibit ‘wet’ or summer dominated flows with high flow constancy and predictability and the variable summer extremely intermittent flow regime rivers, which also display a high degree of flow predictability but the seasonality of flow is much weaker (Pusey et al. 2009).

The majority of rivers in the Kimberley region have not been subject to anthropogenic hydrological alteration. However, the Ord River is the main exception. The construction of the Ord River Dam has changed downstream flows, since its construction in the 1960s. The Ord Dam impoundment, Lake Argyle, has a capacity of 10,700,000 ML and regulates flows from 90% of the Ord catchment (Doupé and Pettit 2002). Flow regulation has significantly reduced mean and peak flows in the lower Ord River, while increasing base flows during the dry season (Start and Handasyde 2002).

Methods

Initially the river networks of the 10 drainage basins in the Kimberley region were characterised according to the procedure outlined by Thoms et al. (2018). A summary of the approach is outlined here. The river network of the Kimberley Region was digitally derived from 1:250,000 scale topographic maps using a Geographic Information System (ArcGIS 9.3). All drainage lines greater than 10 km in length were included in the dataset and these were cross-checked against the 2009 LANDSAT 4–5 TM satellite imagery of the region at a map scale of 1:100,000. A series of sites were then created along the entire river network at approximately 10 km intervals. Each site became the location for the extraction of 15 geomorphic variables that describe the physical character of the riverine landscape. This was done using a suite of semi-automated GIS tools (Thoms et al. 2018). These 15 variables represent data at three spatial scales; region (> 102 km), valley (101 km) or channel (< 101 km) (Table 1); and have been shown to influence the physical character of riverine landscapes (Leopold and Wolman 1957; Schumm 1977; Petts and Amoros 1996).

Table 1 List of variables used in river network characterisation and the spatial scale at which each operates
Full size table

A suite of multivariate statistics was employed to analyse this large data set (3418 sites by 15 variables) to identified groups of sites of similar morphology. First, a cluster analysis was undertaken using the flexible unweighted pair-group method with arithmetic averages (Flexible-UPGMA) fusion strategy. For this the Gower association measure was used because it is a range-standardised measure recommended for physical data (Belbin 1993). The resultant dendrogram grouped sites of similar physical character. The optimum number of groups selected from the dendrogram was determined by examining the relationship between the number of groups and their level of association. The first major inflexion in this relationship was selected as the optimum number of groups as recommended by Quinn and Keough (2002). This statistical grouping of sites, with similar morphological character, equate to the identification of FPZs (Thoms et al. 2018). This self-emerged statistical grouping of sites was then arrayed back onto the river network to produce a morphological characterisation of the river network of the Kimberley region. Second, to further elucidate the grouping of sites into FPZs, a semi-strong-hybrid multidimensional scaling ordination (MDS) was performed on the data. Sites were arrayed in ordination space and then an ANalysis Of SIMilarity (ANOSIM) was performed to assess statistical differences between groups of sites (FPZs). Third, a SIMilarity PERcentage (SIMPER) analysis was undertaken to determine the contribution of the 15 geomorphic variables to the within group similarity: thus, identifying the variables important in creating the observed similarity of each FPZ. Identification of these variables were used to construct a FPZ nomenclature for the Kimberley river network.

Boundaries; those locations within the river network where different FPZs join, were then identified. The relative ‘strength’ of expected boundaries within the river network of the Kimberley region was determined as the distance between the group centroids, or centroid distances, of the FPZs when arrayed in ordination space. Centroid distances between groups in ordination space are a direct measure of the strength of difference between groups (Quinn and Keough 2002). Thus, in the context of this study, longer centroid distances represent increased morphological differences between potentially adjoining FPZs and thus stronger boundary conditions within the river network. Boundary strengths were then classified into one of six ‘strength classes’ according to Table 2.

Table 2 Boundary strength classes determined a priori based on the ordination distance between river types group centroids
Full size table

The spatial distribution of boundary strengths throughout the Kimberley river network was analysed against a set of environmental variables hypothesised to influence the location of boundaries. These were lithology type, proximity to lithological boundaries, elevation, slope, and occurrence at confluences. Each of which has been inferred to influence boundaries or critical transition zones (Benda et al. 2004; Phillips 2008; Rice et al. 2008). First, the proportion of all boundaries in each of the five lithology types identified in the Kimberley was calculated. These proportions were used as the ‘expected’ proportions in subsequent analyses to account for the uneven distribution of lithology types throughout the Kimberley. The proportion of boundaries in each lithology type was then calculated individually for the boundary strength classes; this was the ‘observed’ proportion for each strength class. The observed/expected ratio was then calculated for each lithology type within each strength class. Second, a 10 km buffer around lithological boundaries throughout the Kimberley region was established and the proportion of boundaries in each strength class that occurred within this buffer was calculated. The distance of 10 km was chosen because this was equivalent to the sampling interval along the river network and represented the minimum sampling resolution. Third, spatial patterns in boundary distributions based on topography were investigated. The existence of a relationship between boundary strength and either elevation or slope was investigated using least squares regression in SPSS. Factors that returned a significant result (p < 0.05) with r2 > 0.8 were considered influential on the distribution of boundary strength. Finally, the proportion of boundaries in each strength class that occurred at network confluences was calculated, as well as the proportion of all confluences throughout the Kimberley river network that were found to be boundaries using the approach adopted here.

Results

River characterisation of the Kimberley region

A total of 35,746 kms of river network across the Kimberley region was analysed. Eleven groups of sites emerged as the optimum number of groups from the cluster analysis, explaining 83% of the similarity between sites. These 11 statistical groups were taken to represent 11 distinct FPZs, each having a similar morphological character that differed from one another. This difference between FPZs was confirmed by the ANOSIM (Global R = 0.749, p < 0.001). Moreover, each FPZ had a unique set of geomorphic variables contributing to the within group similarity, as determined from the SIMPER results (Fig. 2). Overall, valley-scale variables were the dominant contributor to within group similarity for all FPZs, albeit with different contributions (Fig. 2). Valley trough width contributed > 40% of the within group similarity for FPZs 7 and 11, which were located in the lower-most regions of the different sub-catchments and these were associated with broad valleys and extensive floodplain surfaces. Whereas the ratio of valley width to valley trough width was the dominant contributor to within group similarity of FPZs 1, 3 and 4 (Fig. 2). Based on the SIMPER results a nomenclature for the Kimberley river characterisation is given in Table 3.

Fig. 2
figure 2

Percentage contribution of each geomorphic variable to the within group similarity of each Functional Process Zone – FPZ - (from SIMPER). See Table 1 for description of variables. Group 9 contained only one segment in the river network, thus SIMPER analysis was not possible

Full size image
Table 3 Description, abbreviation and total stream length of each of the 11 Functional Process Zones identified throughout the Kimberley river network
Full size table

Overall, five FPZs were abundant in the Kimberley river network; the Headwater zone (Hw), Upland Moderate Slopes zone (UpLMnd), Midland Moderate Slopes zone (MidMnd), Mid to Lowland Gorges zone (GMidLow), Mid to Lowland Anabranching zone (MidLowAnb), and the Single Channel Broad Valley Lowland zone (LowMnd), which contributed to a combined stream length of over 30,000 km or 84% of the total river network (Fig. 3). The headwater zone (Hw) was the most abundant FPZ, in terms of total stream length, constituting 9701 km or 27% of the entire river network of the Kimberley (Fig. 3). Two FPZs, the Sinuous Gorge zone (GHSin) and the Broad Valley Constrained Trough zone (BrdValNrwTr) were rare, with a length of 12 and 72 km of river network, respectively (Fig. 3).

Fig. 3
figure 3

Spatial distribution of the 11 Functional Process Zones (FPZ) identified throughout the Kimberley river network

Full size image

Marked spatial patterns in the distribution of FPZs were evident across the Kimberley region (Fig. 3). FPZs Hw and UpLMnd were found predominantly in the upper sections of most drainage basins, with the former being widespread throughout the Kimberley Plateau, while the latter was more abundant in the eastern Kimberley but less common in other areas and completely absent in the far northern and south-western parts of the region (Fig. 3). By comparison, the MidLowAnb, LowMnd, and BrdLowAnb FPZs were strongly associated with the lower sections of most rivers, particularly in the Fitzroy, Lennard and lower Ord basins (Fig. 3). Moreover, some FPZs were uniquely associated with particular physiographic areas of the Kimberley region. FPZs MltChanMid and GUp, for example, were relatively common in the King Leopold and Durack Ranges (Figs. 1 and 3) but uncommon elsewhere.

River network boundaries

A total of 1410 boundaries were identified throughout the Kimberley river network. This represents an average of one morphological boundary or potential discontinuity for every 25 km of stream length within this riverine landscape (Fig. 4). By comparison, there were 914 tributary junctions. Of the possible 55 boundary types among the 11 FPZs, only 41 of these were observed in the Kimberley river network. The strength of the identified boundaries ranged from 0.0789 to 0.5425 (Table 2) with the most frequent boundary occurring between the Hw and UpLMnd FPZs. Statistically this was also the weakest boundary with a strength of 0.0789. The strongest boundary was between the GUp and BrdLowAnb FPZs, with a strength of 0.5452, and this occurred only once in the region, in the upper Ord basin (Fig. 4D). The distribution of boundary strengths was positively skewed, with boundary classes 2 and 3 being dominant (Fig. 4).

Fig. 4
figure 4

The spatial distribution of boundaries throughout the Kimberley river network and the five broad lithology types. Boundaries are displayed by strength class: A class 1, B class 2, C class 3, and D classes 4 to 6

Full size image

Spatial distribution and environmental drivers of river network boundaries

Distinct patterns in the distribution of the different boundary strengths occurred throughout the Kimberley river network. Boundaries in class 1 were mainly located in the central and eastern parts of the Kimberley (Fig. 4), while boundaries in classes 2 and 3 had a relatively uniform distribution over most of the region (Fig. 4B and C), with the exception of the central plateau, where those in class 3 were scarce (Fig. 4B). By comparison, stronger boundaries in classes 4, 5 and 6 were relatively rare; representing only 2.2% (n = 31) of the total number of boundaries identified.

Boundaries in the Kimberley river network mostly occurred in the sedimentary lithology class (65.4%), which is unsurprising given this is the region’s dominant lithology. However, the strength classes of boundaries were disproportionately associated with other lithology types. Boundaries in strength class 1 were associated more with the less common lithology types than with sedimentary rocks (Fig. 4A). The proportion of boundaries in class 1 that were associated with the granulite-facies metamorphics lithology type was more almost four times that which would be expected based on the distribution of all nodes among different lithology types (Fig. 5). Stronger boundaries (classes 4, 5 and 6) were underrepresented in areas of sedimentary rocks and occurred mainly in other lithology types (Figs. 4D and 5). Over 25% of the boundaries in class 4 occurred on mafic-ultramafic intrusives, dolerites and gabbros, compared to only 6% of all boundaries in the region. The strongest boundary, of which there was only one in class 6, occurred within the mafic and felsic volcanic lithology type, which contained only 20% of all boundaries in the river network, representing a five-fold increase in the observed/expected ratio. The number of boundaries in class 5 that occurred on mafic and felsic volcanic lithology class was also more than twice the expected (Fig. 5).

Fig. 5
figure 5

The observed/expected ratio of each boundary strength class within each lithology type

Full size image

73% of all boundaries that occurred throughout the Kimberley river network were within 10 km of a lithological boundary (cf. Table 4). This proportion was accentuated for boundary classes 1 and 4, in which 93 and 96% of boundaries occurred within 10 km of a lithological boundary, respectively (Table 4). Boundaries in the two strongest classes (5 and 6) did not occur within 10 km of a lithological boundary. Associations between boundary strength and topography varied. A significant increase in boundary strength with decreasing elevation was observed (F = 46.320; d.f. = 1, 1409; p < 0.000); however, the relationship was weak (r2 = 0.032). No relationship between boundary strength and slope throughout the Kimberley region was evident (F = 1.469; d.f. = 1, 1409; p = 0.226).

Table 4 Total number and proportions of boundaries in each strength class that occurred near lithological boundaries and at river network confluences
Full size table

Of the 1410 boundaries identified between the 11 different FPZs in the Kimberley river network, 32% occurred at confluences (Table 4), while the remaining 68% were located elsewhere along the river network. However, the majority of strong boundaries (classes 5 and 6) occurred at confluences (Table 4). A total of 914 confluences exist throughout the Kimberley river network used in this study. Thus, half of the total number of confluences in the region were not found to be boundaries between morphologically distinct FPZs using the characterisation approach adopted here (Fig. 6).

Discussion

Our analysis shows that the river network of the Kimberley region is discontinuous, as indicated by our finding that a geomorphic boundary occurs every 25 km of river, on average. Boundaries or discontinuities within the river network result from a range of lithological, hydrological, and geomorphological processes. For example, the strongest boundary we observed was where an upland gorge FPZ transitioned into a lowland anabranching channel FPZ. This boundary manifest as a transition from a highly constrained channel setting with high flow velocities to one of multiple unconstrained channels with lower velocities and a finer river-bed substratum. Other weaker but far more numerous boundaries were observed, for example, where upland moderate slope FPZs transitioned into mid to lowland anabranching channels FPZ—representing a decrease in slope and a transition from single to multi-channelled river sections. The presence of geomorphic boundaries shapes the discontinuous character of riverine landscapes in the Kimberley, which may support the high freshwater biodiversity of the region (Pueschel 2019).

The presence of geomorphic boundaries influences the longitudinal continuum/discontinuum within river networks (Naiman et al. 1988; Phillips 2008). We suggest that the degree of connectivity or character of the river continuum, is determined by the contrast in geomorphic character between two adjacent FPZs. A range of connectivities, as indicated by the magnitude and frequency of boundary strengths, were calculated for the Kimberley river network. Most of the boundaries in the Kimberley (51%) were weak, i.e., had a boundary strength of class 1 or 2. By contrast, less than 1% of the boundaries in the Kimberley had strengths of five or greater. The distribution of boundary strengths, as characterised by their magnitude and frequency, determines the overall longitudinal continuum in river networks. The nature of this distribution and the spatial organisation of geomorphic boundaries will dictate fragmentation – continuum/discontinuum of riverine landscapes.

Our analysis reveals multiple associations between environmental drivers, the spatial organisation and, strength of boundaries. Lithological boundaries have a strong influence on the location of river network boundaries in the Kimberley region, particularly those between relatively similar FPZs (i.e., ‘weaker’ boundaries). In contrast, the strongest boundaries are mostly associated with tributary confluences and are disproportionately present in two uncommon lithology types: (1) mafic and felsic volcanics, and (2) mafic-ultramafic intrusives, dolerites and gabbros. These strong boundaries may be caused by greater physical heterogeneity in these lithologies because of their rock properties. While tributary confluences appear to create the strongest river network boundaries, they do not determine the number and location of most boundaries. Similarly, only half of all tributary confluences were found to be boundaries using our method, suggesting that a broader perspective on river network discontinuities is necessary. Our results also indicate that elevation and slope are not main drivers of the location or strength of river network boundaries in the Kimberley region.

Scale is an important factor affecting our results and assessments of river network discontinuities in general. The NDH of Benda et al. (2004) focuses on tributary confluences and in particular the ratio of tributary discharge to main stem channel discharge. In this way, the contrast between reaches and the strength of the boundary are determined locally at each confluence. In contrast, geomorphic boundaries and their strengths are determined by a multi-basin network-scale characterisation of FPZs. Our results reveal that, at this larger scale, confluences become less important for determining boundaries. Another effect of scale is the scale of variables used to determine the presence and strength of boundaries. In the NDH, this is largely based on within-channel variables associated with discharge (Benda et al. 2004). Our approach is based on regional, valley, and channel planform variables, so we cannot determine any within-channel effects of the boundaries. The sampling resolution of our approach (one site approximately every 10 km along the network) is appropriate for a multi-basin analysis across tens of thousands of stream kilometres, but a finer sampling resolution is required for more detailed analyses. This is possible with the method employed by simply creating more frequent sampling sites in GIS. The types of boundaries that can be observed also depends on scale and the number of groups taken from the hierarchical cluster analysis—results ultimately depend on the contrast between two distinct categorical groups. An alternative ‘continuous’ as opposed to categorical approach would be to quantify the multidimensional distance (using the same 15 variables or others) between sites consecutively downstream, rather than grouping all reaches into FPZs first.

While longitudinal continuity is an essential feature of river networks (Ward and Stanford 1995; Cote et al. 2009), discontinuity is also (Poole 2002). Boundaries create modularity, which protects systems against contagious or catastrophic disturbance (Ash and Newth 2007), and the maintenance of boundaries in river networks provides heterogeneity and refugia (Grant et al. 2007). Connections among patches (e.g., functional process zones and stream reaches) of different morphology helps maintain natural discontinuity and diversity, which likely bolster the resilience of river networks. Understanding the locations of boundaries and what causes them throughout river networks allows us to explore questions such as: where in the network are geomorphic processes disrupted? And where might greater physical heterogeneity occur? Similarly, knowing the relative strengths of boundaries might indicate the likelihood of disruptions in downstream processes or hotspots of physical or biological diversity (e.g., Benda et al. 2004; Czuba and Foufoula-Georgiou, 2015). Most focus so far in answering these questions has been on tributary confluences (cf. Benda et al. 2004; Rice et al. 2008), and research suggests that the probability of morphological and biological effects of confluences can be predicted by the ratio of tributary discharge to main stem channel discharge (Benda et al. 2004; Kiffney et al. 2006). Studies have shown that the richness of fish species throughout a river network are weakly correlated with confluence location (Hitt and Angermeier 2008; Boddy et al. 2019). This focus on tributaries disregards potentially important boundaries at other locations throughout river networks, and our results show that there can be many more. Thus, novel quantitative approaches such as ours provide a more complete picture of discontinuity in river networks (Fig. 6).

Fig. 6
figure 6

Comparison of our quantitative river characterisation approach and a tributary confluence-based approach to determining boundaries in river networks. The number of boundaries is 1410 under our approach and 914 based on confluences. Our approach reveals additional boundaries not at confluences (a) and omits around half of the confluences in the network (c)

Full size image