Optimal Channel Networks accurately model ecologically-relevant geomorphological features of branching river networks
Drainage area and branching ratio: a matter of scaleGeomorphological and ecological viewpoints on river networks generally differ owing to discordant definitions of the fundamental unit (the node) used to analyze them. From a geomorphological perspective, the determination of a river network entails the definition of an observational scale. Real river networks can be extracted from digital elevation models (DEMs) via algorithms for flow direction determination such as D8 (i.e., each pixel drains towards the lowest of its 8 nearest neighbors53). After the outlet location has been specified (and hence the upstream area A spanned by the river network), the first observational scale required is thus the pixel length l of the DEM, which defines the extent of a network node. A second scale is then needed to distinguish the portion of the drainage network effectively belonging to the channel network. The simplest but still widely used method53 defines channels as those pixels whose drainage area exceeds a threshold value AT. Hydrologically based criteria to determine the appropriate value for AT exist54; however, for the sake of simplicity, we here consider AT as a free parameter.BBTs and RBNs are random constructs, and as such they do not satisfy the optimality criterion of minimizing total energy expenditure, which is the fundamental physical process shaping fluvial landscapes. Furthermore, neither of these networks is a spanning tree, which is a key attribute of real fluvial landforms10: in fact, in both BBTs and RBNs, the extent of the drained domain is not defined. As a result, the drainage area at an arbitrary network node cannot in principle be attributed, unless by using the number of upstream nodes as a proxy. This has practical implications from an ecological viewpoint because drainage area is the master variable controlling several attributes of a river, such as width, depth, discharge, or slope3,55, which in turn impact habitat characteristics and the ecology of organisms therein56.In BBTs and RBNs, branching probability p has been defined35,38,45,46,47 as the probability that a network node is branching, i.e. connected to two upstream nodes. As such, the branching probability of a realized river network (be it a real river or a synthetic construct) could be evaluated as the ratio between the number of links NL constituting a network and the total number of network nodes N; if a unit distance between two adjacent nodes is assumed, the denominator equals the total network length. We note that the former definition of branching probability only holds in the context of the generation of a synthetic random network; it is in fact improper to refer to a “probability” when analyzing the properties of a realized river network. We clarify this aspect by introducing the concept of branching ratio pr for the latter definition (pr = NL/N). Moreover, in the case of BBTs, p and pr do not coincide (see Methods). Importantly, p and pr have no parallel in the literature on fluvial forms, nor do they refer to any of the well-studied measures of rivers’ fractal character.The choice of different observational scales for the same drainage network results in different values of NL and N, and hence of pr. Remarkably, the very same drainage network can result in river networks that virtually assume any value of pr (ranging from 0 to 1) and N (up to the upper bound A) depending on the choice of AT and A (the latter corresponding to a given l value when measured in the number of pixels; Fig. 1d–i); networks with low AT/A ratios result in high N (Fig. 2a), while networks with low AT result in high pr (Fig. 2b). Furthermore, pr does not identify the inherent (i.e., scale-independent) branching character of a given river network in relation to other river networks. In fact, by extracting different river networks at various scales (i.e., various AT values) and assessing the rivers’ rank in terms of pr, one observes that rivers that look more “branching” (i.e., have higher pr) than others for a given AT value can become less “branching” for a different AT value (Fig. 3). We therefore conclude that branching probability is a non-descriptive property of a river network, which by no means describes its inherent branching character, and depends on the observational scale.Fig. 2: Variation of N and pr as a function of observational scales for OCNs and real river networks.a Expected value of number of network nodes N as a function of threshold area AT and total drained area A (from Eq. (1)); the white dots indicate the values of AT and A used to generate the OCNs used in this analysis. b Expected value of branching ratio pr as a function of AT and A (from Eq. (1)); symbols as in a.Full size imageFig. 3: Values of branching ratio as a function of AT for the 50 real river networks analyzed in this study.a Natural values of pr in logarithmic scale. b z-normalized branching ratios (i.e., for each AT value, values of pr are normalized so that they have null mean and unit standard deviation), which better shows how rivers rank differently in terms of pr for different observation scales (i.e., AT). Lines connect dots relative to the same river. For visual purposes, rivers that rank first, second, second-to-last or last in at least one of the AT groups are displayed in colors; the other rivers are displayed in grey.Full size imageScaling is also crucial when looking at river networks from an ecological perspective. In this case, the relevant scale determining the dimension l of a node is the extent of habitat within which individuals (due to e.g. physical constrains) can be assigned to a single population57,58; the riverine connectivity and ensuing dispersal among these populations give rise to a metapopulation at the river network level. The specific spatial scale largely depends on the targeted species (e.g. being larger for fish than for aquatic insects), and it is conceivably much larger than (or, at least, it has no reason to be equal to) the pixel size of the DEM on which the river network is extracted. Since the evaluation of pr depends on the number of nodes N, which, in turn, is defined based on the scale length l, the resulting pr of a river network under this perspective would depend on the characteristics of the target taxa, which is inconsistent with the alleged role of pr as a scale-invariant property of river networks.Note also that using the ecological definition of l (i.e., spatial range of a local population) to discretize a real river network into N nodes, and from there calculate the branching ratio pr = NL/N, is problematic. Indeed, this would imply an elongation of all links shorter than l (which constitute a non-negligible fraction of the total links, under the assumption of exponential distribution of link lengths51), hence preventing a correct estimation of the connectivity patterns (i.e., distances between nodes) and the resulting ecological metrics of the river network (see section Ecological implications).From an ecological perspective, it could be reasonable to consider AT as a parameter expressing how a particular taxon perceives the suitable landscape, rather than a value to be determined from geomorphological arguments: for instance, large fishes inhabit wide and deep river reaches, and do not access small headwaters56. In this case, imposing a large AT would result in a coarser, less branching network constituted by few main channels (Fig. 1f, i), which could mimic the potentially available habitat for such species. Conversely, aquatic insects inhabit also small headwaters17,59, therefore their perceived landscape would resemble the finely resolved networks of Fig. 1d, g, characterized by low AT and higher (apparent) pr.Topology and scaling of river networks and random analoguesTo verify the topological (i.e., Horton’s laws on bifurcation and length ratios) and scaling (i.e., probability distribution of drainage areas) relationships of the different network types, we extracted from DEMs 50 real river networks encompassing a wide range of drainage areas (Fig. 4), and we generated 50 OCNs, 50 RBNs and 50 BBTs of comparable size (see Methods).Fig. 4: Location of real river basins used in the analysis.River basins are shown in dark grey; countries in light grey. Rivers’ numbering is sorted in ascending order according to drainage area values.Full size imageTypical values3,7,60 for the bifurcation ratio RB lie between 3 and 5, while length ratios (RL) range between 1.5 and 3.5. As expected, we observed that the real rivers and OCNs used in our analysis have RB and RL values within the aforementioned ranges (Fig. 5a, b). The same is true for RBNs, while the RB and RL values found for BBTs are lower than the typical ranges. This finding holds regardless of the scale (subsumed by AT) at which real river networks and OCNs are extracted (Supplementary Figs. 1 and 2). Remarkably, BBTs fail to satisfy Horton’s laws despite the statistical inevitability of such laws for any network argued by ref. 61. To this regard, we note that the networks analyzed by ref. 61 did not include constructs where all paths from the source nodes to the outlet have the same length, which is the defining feature of BBTs (Fig. 1a).Fig. 5: Comparison of topological and scaling properties of the different networks.a Scaling of number of network links Nω as a function of stream order ω for the various network types (rivers and OCNs obtained with AT = 20 pixels; RBNs and BBTs derived accordingly – see Methods). b Mean link length Lω (in units of l) as a function of ω. Networks used are as in panel a. c Scaling of drainage areas: probability P[A ≥ a] to randomly sample a node with drainage area A ≥ a as a function of a. The displayed trend lines are fitted on the ensemble values for the 50 network replicates, by excluding nodes with drainage area larger than 2000 pixels (cutoff value marked with a black solid line). The scaling coefficients β reported correspond to the slopes of the fitted trend lines. Extended details on all panels are provided in the Supplementary Methods.Full size imageWhile the power-law scaling of areas in OCNs (Fig. 5c) has an exponent β ≈ 0.45 that closely resembles the one found for the real rivers (β ≈ 0.46) and within the typically observed range8,10 β = 0.43 ± 0.02, drainage areas of RBNs scale as a power law with an exponent β ≈ 0.51, which departs from the observed range. Conversely, BBTs do not show any power-law scaling of areas. Scaling exponents of drainage areas fitted separately for each real river network yielded values in the range 0.36÷0.57 (Supplementary Table 1). In particular, we observed that these values tend to the expected range β = 0.43 ± 0.02 for increasing values of A, expressed in number of pixels (Supplementary Fig. 3), hence implying that highly resolved catchments are required in order to properly estimate β. Interestingly, the observed values of Horton ratios and scaling exponent β for RBNs are compatible with the values RB = 4, RL = 2, β = 0.5 predicted for Shreve’s random topology model3,60,62, which is actually equivalent to a RBN with infinite links.Ecological implicationsWe compared the different network types via two metrics that express the ecological value of a landscape for a metapopulation: the coefficient of variation of a metapopulation CVM and the metapopulation capacity λM. The coefficient of variation of a metapopulation63 is a measure of metapopulation stability (a metapopulation being more stable the lower CVM is), while the metapopulation capacity42,64 expresses the potential for a metapopulation to persist in the long run (persistence being more likely the higher λM is). Both measures are among the most universal metrics describing dynamics of spatially fragmented populations24,40. In order to assess the impact of the two landscape features mostly affecting metapopulation dynamics, i.e. spatial connectivity and spatial distribution of habitat patches, we calculated these metrics for the four network types under two different scenarios: uniform (CVM,U, λM,U) and non-uniform (CVM,H, λM,H) spatial distribution of habitat patch sizes. In the first scenario, CVM,U and λM,U assess stability and persistence (respectively) of a metapopulation solely based on pairwise distances between network nodes; in the second scenario, CVM,H and λM,H depend on the interplay between pairwise distances and spatially heterogeneous habitat availability (namely, downstream nodes being larger than upstream ones).We found that the values of CVM (be it derived with uniform (CVM,U) or nonuniform (CVM,H) distributions of patch sizes) obtained for OCNs match strikingly well those of real rivers (Fig. 6). These CVM values are consistently lower than those found for RBNs, while values of CVM for BBTs are even higher. Notably, this result holds for different values of AT (and hence different pr values) at which real rivers and OCNs are extracted (Fig. 6a–c; g–i), and for values of mean dispersal distance α (see Methods) spanning multiple orders of magnitude (Supplementary Figs. 4–7).Fig. 6: Comparison of values of metapopulation metrics across river network types and observational scales (AT).a–c CVM,U. d–f λM,U. g–i CVM,H. j–l λM,H. Boxplot elements are as follows: center line, median; notches, (pm 1.58cdot {{{{{{{rm{IQR}}}}}}}}/sqrt{50}), where IQR is the interquartile range; box limits, upper and lower quartiles; whiskers, extending up to the most extreme data points that are within ±1.5 ⋅ IQR; circles, outliers. Metapopulation metric values were obtained by setting α = 100 l. Note that in Eq. (1), given A = 40, 000, AT = 20 results in E[N] ≈ 4574, E[pr] ≈ 0.228; AT = 100 yields E[N] ≈ 2231, E[pr] ≈ 0.098; AT = 500 results in E[N] ≈ 1088, E[pr] ≈ 0.042.Full size imageFor a constant α value, the CVM of real rivers, OCNs and RBNs decreases as the resolution at which the network is extracted increases (i.e., AT decreases; see Fig. 6 and Supplementary Figs. 4–7). This is expected63, since a decrease in AT corresponds to an increase in N (Fig. 2a), leading to a decrease in CVM. Indeed, a larger ecosystem, constituted of more patches, has the potential to include a larger (and more diverse) number of subpopulations, which increases stability at a metapopulation level through statistical averaging–a phenomenon widely known as the portfolio effect65. We also found that BBT networks do not generally follow the above-described pattern of decreasing CVM with increasing N; rather, the CVM of BBTs increases with N when the mean dispersal distance α is set to intermediate to high values (Fig. 6 and Supplementary Figs. 5–7), and only when α is very low (e.g. α = 10 l as in Supplementary Fig. 4) and a uniform patch-size distribution is assumed does CVM,U follow the expected decreasing trend with increasing N.However, we need to warn against the conclusion that river networks with higher values of pr (and hence lower AT, see Fig. 2b) are inherently associated with higher metapopulation stability. Indeed, our result was obtained by changing the scale at which we observed the same river networks, and not by increasing the river networks’ size. If the number of network nodes (and, consequently, the branching ratio pr) is determined by the scale at which the landscape is observed, one cannot directly assume that any of such nodes is a node (or patch) in the ecological sense, i.e. the geographical span of a local population: the extent of such patches should be determined based on the mobility characteristics of the focus species, and should be independent of the scale at which the river network is observed. In contrast, we note that, if different river networks spanning different catchment areas (say, in km2) are compared, all of them extracted from the same DEM (same l and same AT in km2), then the larger river network will appear more branching (i.e., have larger pr). Indeed, by selecting catchments with larger A (in km2) for fixed l and AT (in km2), one moves towards the top-left corner of Fig. 2a, b (i.e., perpendicular to the level curves AT/A). The apparent higher “branchiness” of the river network with larger A will result in lower values of CVM; however, the higher metapopulation stability of the larger network will not be due to its (alleged) inherent more branching character, but only dictated by its larger habitat availability.We observed that metapopulation capacity λM values of OCNs (be it evaluated under uniform (λM,U) or non-uniform (λM,H) patch-size distribution assumption) are the closest to those of real rivers, while RBNs (and even more so BBTs) generally overestimate λM with respect to real rivers and OCNs (Fig. 6d–f; j–l). This result holds irrespective of the choice of AT and for intermediate to high values of α (Supplementary Figs. 5–7). When the mean dispersal distance is instead set to very low values (α = 10 l – Supplementary Fig. 4) and the river network is extracted at a high resolution (i.e., low AT), the metapopulation capacity of OCNs under assumption of uniform patch-size distribution (λM,U) is underestimated with respect to that of real rivers. A likely explanation for this apparent mismatch is that, for low values of AT, the number of nodes N tends to be somewhat higher for the extracted river networks used in this analysis than for OCNs (Supplementary Fig. 8), and the effect of the different dimensionality of real rivers and OCNs in the metapopulation capacity estimation tends to be more evident as the mean dispersal distance decreases. Interestingly, such mismatch is absent when a non-uniform patch size distribution is assumed, as λM,H values for OCNs match those for real rivers regardless of the mean dispersal distance value and the river network resolution (Fig. 6; Supplementary Figs. 4–7).The OCN construct encapsulates both random and deterministic processes, the former related to the stochastic nature of the OCN generation algorithm, and the latter pertaining to the minimization of total energy expenditure that characterizes OCN configurations. As such, OCNs reproduce the aggregation patterns of real river networks. From an ecological viewpoint, this implies that both pairwise distances between nodes and the distribution of patch sizes (expressed as a function of drainage areas, or of a proxy thereof such as the number of nodes upstream) are much closer to those of real networks than is the case for fully random synthetic networks as BBTs and RBNs. In particular, BBTs and (to a lesser extent) RBNs tend to underestimate pairwise distances with respect to real rivers and OCNs, as documented by a comparison of mean pairwise distances across network types (Supplementary Fig. 9a–c). Our analysis shows that the connectivity structure of these random networks (subsumed by the matrix of pairwise distances) is too compact with respect to that of real rivers, which leads to an overestimation of the role of dispersal in increasing the ability of a metapopulation to persist in the long run, but also an increased likelihood of synchrony among the different local populations, which results in higher instability.Comparison of patch size distributions among the network types expressed in terms of CVM,0 (i.e., the portion of CVM,H that uniquely depends on the distribution of patch sizes and not on pairwise distances) shows that, while for coarsely resolved networks (AT = 500) no clear differences in CVM,0 emerged, for highly resolved networks (AT = 20) BBTs heavily underestimate the CVM,0 of real rivers and OCNs, while RBNs slightly overestimate it (Supplementary Fig. 9d–f). As a result of the interplay of differences in distance matrices and patch size distributions, BBTs and (to a lesser extent) RBNs generally tend to overestimate the coefficient of variation of a metapopulation and the metapopulation capacity of real rivers and OCNs in both scenarios of uniform and non-uniform patch size distribution. The only exception to this trend occurs for the metapopulation capacity λM,H of very large BBTs (corresponding to AT = 20) in the case of very high dispersal distances (α = 1000 l – Supplementary Fig. 7): here, the patch-size effect (i.e., underestimation of CVM,0) predominates over the distance effect (i.e., overestimation of mean dij), resulting in an underestimation of λM,H with respect to real rivers and OCNs.Our results were derived under a number of simplifying assumptions. In particular, we acknowledge that, while the distance matrix of a landscape and the distribution of patch sizes have in general important implications for metapopulation dynamics, other factors not considered here, such as Euclidean between-patch distance48, fat-tailed dispersal kernel66 and density-dependent dispersal67 could also play a relevant role in this respect. However, it needs to be noted that, especially with regards to the assessment of the Moran effect in metapopulation synchrony (i.e., increased synchrony in local fluvial populations that are geographically close but not flow-connected48), the use of OCNs allows integration of Euclidean distances in a metapopulation model, while this is not possible for RBNs and BBTs, where Euclidean distances are not defined. Moreover, if a larger degree of realism is required for a specific ecological modelling study, such as heterogeneity in abiotic factors (e.g. water temperature or flow rates), the use of OCNs as model landscapes allows a direct integration of these variables, as they can conveniently be expressed as functions of drainage area3,55. In contrast, this is not possible for RBNs or BBTs, because only OCNs verify the scaling of areas (Fig. 5c), while RBNs and BBTs lack a proper definition of drainage areas.Our comparison of synthetic and real river networks showed that riverine metapopulations are more stable and less invasible than what would be predicted by random network analogues. Conversely, the use of OCNs as model landscapes allows capturing not only the scaling features of real rivers, but also drawing ecological conclusions that are in line with those that could be observed in real river networks. We thus support the use of OCNs as analogues of real river networks in theoretical and applied ecological modelling studies. While we found that BBTs are highly inaccurate in reproducing ecological metrics of real river networks and should be therefore discarded altogether in future modelling applications, RBNs show a certain degree of similarity with OCNs and real river networks in this respect; moreover, RBNs (as is the case for any random tree61) satisfy Horton’s laws on bifurcation and length ratios. A relevant advantage of RBNs over OCNs is that their generation algorithm is at least one order of magnitude faster49. Therefore, we acknowledge that RBNs could be considered as a suitable surrogate for real river networks as null models in cases where a large number of network replicates is required. To this end, we encourage researchers exploiting synthetic river networks (whether they be OCNs or RBNs) to always clarify the observational scales (that is, total area drained, size of a node, area drained by a headwater) subsumed by the synthetic network and which give rise to a certain complexity measure (i.e., branching ratio). Only in such a way could the predictions from these studies be compared with real river networks.In conclusion, our results advocate a tighter integration between physical (geomorphology, hydrology) and biological (ecology) disciplines in the study of freshwater ecosystems, and particularly in the perspective of a mechanistic understanding of drivers of persistence and loss of biodiversity. More
