Optimal Channel Networks accurately model ecologically-relevant geomorphological features of branching river networks

Drainage area and branching ratio: a matter of scale

Geomorphological 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 neighbors⁵³). 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 method⁵³ defines channels as those pixels whose drainage area exceeds a threshold value A_T. Hydrologically based criteria to determine the appropriate value for A_T exist⁵⁴; however, for the sake of simplicity, we here consider A_T 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 landforms¹⁰: 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 slope^3,55, which in turn impact habitat characteristics and the ecology of organisms therein⁵⁶.

In BBTs and RBNs, branching probability p has been defined^{35,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 N_L 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 p_r for the latter definition (p_r = N_L/N). Moreover, in the case of BBTs, p and p_r do not coincide (see Methods). Importantly, p and p_r 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 N_L and N, and hence of p_r. Remarkably, the very same drainage network can result in river networks that virtually assume any value of p_r (ranging from 0 to 1) and N (up to the upper bound A) depending on the choice of A_T and A (the latter corresponding to a given l value when measured in the number of pixels; Fig. 1d–i); networks with low A_T/A ratios result in high N (Fig. 2a), while networks with low A_T result in high p_r (Fig. 2b). Furthermore, p_r 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 A_T values) and assessing the rivers’ rank in terms of p_r, one observes that rivers that look more “branching” (i.e., have higher p_r) than others for a given A_T value can become less “branching” for a different A_T 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.

Scaling 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 population^57,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 p_r depends on the number of nodes N, which, in turn, is defined based on the scale length l, the resulting p_r of a river network under this perspective would depend on the characteristics of the target taxa, which is inconsistent with the alleged role of p_r 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 p_r = N_L/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 lengths⁵¹), 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 A_T 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 headwaters⁵⁶. In this case, imposing a large A_T 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 headwaters^17,59, therefore their perceived landscape would resemble the finely resolved networks of Fig. 1d, g, characterized by low A_T and higher (apparent) p_r.

Topology and scaling of river networks and random analogues

To 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).

Typical values^3,7,60 for the bifurcation ratio R_B lie between 3 and 5, while length ratios (R_L) range between 1.5 and 3.5. As expected, we observed that the real rivers and OCNs used in our analysis have R_B and R_L values within the aforementioned ranges (Fig. 5a, b). The same is true for RBNs, while the R_B and R_L values found for BBTs are lower than the typical ranges. This finding holds regardless of the scale (subsumed by A_T) 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. ⁶¹. To this regard, we note that the networks analyzed by ref. ⁶¹ 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).

While 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 range^8,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 R_B = 4, R_L = 2, β = 0.5 predicted for Shreve’s random topology model^3,60,62, which is actually equivalent to a RBN with infinite links.

Ecological implications

We 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 CV_M and the metapopulation capacity λ_M. The coefficient of variation of a metapopulation⁶³ is a measure of metapopulation stability (a metapopulation being more stable the lower CV_M is), while the metapopulation capacity^42,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 populations^24,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 (CV_M,U, λ_M,U) and non-uniform (CV_M,H, λ_M,H) spatial distribution of habitat patch sizes. In the first scenario, CV_M,U and λ_M,U assess stability and persistence (respectively) of a metapopulation solely based on pairwise distances between network nodes; in the second scenario, CV_M,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 CV_M (be it derived with uniform (CV_M,U) or nonuniform (CV_M,H) distributions of patch sizes) obtained for OCNs match strikingly well those of real rivers (Fig. 6). These CV_M values are consistently lower than those found for RBNs, while values of CV_M for BBTs are even higher. Notably, this result holds for different values of A_T (and hence different p_r 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).

For a constant α value, the CV_M of real rivers, OCNs and RBNs decreases as the resolution at which the network is extracted increases (i.e., A_T decreases; see Fig. 6 and Supplementary Figs. 4–7). This is expected⁶³, since a decrease in A_T corresponds to an increase in N (Fig. 2a), leading to a decrease in CV_M. 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 effect⁶⁵. We also found that BBT networks do not generally follow the above-described pattern of decreasing CV_M with increasing N; rather, the CV_M 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 CV_M,U follow the expected decreasing trend with increasing N.

However, we need to warn against the conclusion that river networks with higher values of p_r (and hence lower A_T, 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 p_r) 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 km²) are compared, all of them extracted from the same DEM (same l and same A_T in km²), then the larger river network will appear more branching (i.e., have larger p_r). Indeed, by selecting catchments with larger A (in km²) for fixed l and A_T (in km²), one moves towards the top-left corner of Fig. 2a, b (i.e., perpendicular to the level curves A_T/A). The apparent higher “branchiness” of the river network with larger A will result in lower values of CV_M; 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 A_T 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 A_T), 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 A_T, 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 CV_M,0 (i.e., the portion of CV_M,H that uniquely depends on the distribution of patch sizes and not on pairwise distances) shows that, while for coarsely resolved networks (A_T = 500) no clear differences in CV_M,0 emerged, for highly resolved networks (A_T = 20) BBTs heavily underestimate the CV_M,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 A_T = 20) in the case of very high dispersal distances (α = 1000 l – Supplementary Fig. 7): here, the patch-size effect (i.e., underestimation of CV_M,0) predominates over the distance effect (i.e., overestimation of mean d_ij), 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 distance⁴⁸, fat-tailed dispersal kernel⁶⁶ and density-dependent dispersal⁶⁷ 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-connected⁴⁸), 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 area^3,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 tree⁶¹) 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 faster⁴⁹. 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.

Source: Ecology - nature.com