Channel concavity controls planform complexity of branching drainage networks

The planform geometry of branching drainage networks controls the topography of landscapes and their geomorphic, hydrologic, and ecologic functionality. The complexity of networks' geometry shows significant variability, from simple, straight channels that flow along the regional topographic gradient to intricate, tortuous flow patterns. This variability in complexity presents an enigma, as models show that it emerges independently of any heterogeneity in the environmental conditions. We propose to quantify networks' complexity based on the distribution of lengthwise asymmetry between paired flow pathways that diverge from a divide and rejoin at a junction. Using the lengthwise asymmetry definition, we show that the channel concavity index, describing downstream changes in channel slope, has a primary control on the planform complexity of natural drainage networks. An analytic model and optimal channel network simulations employing an energy minimization principle reveal that landscapes with low concavity channels attain planform stability only with simple network geometry. In contrast, landscapes with high concavity channels can achieve planform stability with various configurations, displaying different degrees of network complexity, including extremely complex geometries. Consequently, landscapes with high concavity index channels can preserve the legacy of former environmental conditions, whereas landscapes with low concavity index channels reorganize in response to environmental changes, erasing the former conditions. Consistent with previous findings showing that channel concavity correlates with climate aridity, we find a significant empirical correlation between aridity and network complexity, suggesting a climatic signature embedded in the large-scale planform geometry of landscapes.

1 Introduction

The planform structure of branching fluvial drainage networks has far-reaching implications for the geomorphic, hydrologic, and ecologic functionality of landscapes (Horton, 1945; Sharp and Malin, 1975; Perron et al., 2008; Willett et al., 2018; Pelletier et al., 2018; Stokes and Perron, 2020; Freund et al., 2023; Liu et al., 2024). This structure, which can be expressed based on its geometric and topological attributes, exhibits significant variation across different regions. In some cases, networks exhibit simple flow paths (Fig. 1a) that generally follow the regional topographic gradient. These flow paths define main drainage basins, draining the main water divide to the mountain front, that are overall similar in shape and size and have a symmetric basin shape with respect to their main trunk (sensu Ramsey et al., 2007). Other networks appear more intricate. These complex networks display tortuous flow paths, asymmetric basin shapes, and varying sizes and shapes of the main basins (Fig. 1b).

Figure 1 Variability of complexity. Examples and schematic representations of paired flow pathways (light blue) originating from two channel heads (yellow squares), diverging from a common divide, and merging at a downstream junction (red square). Hillshade topographies with drainage networks colored in white are displayed in panels (a) and (b). Panel (a) shows an example from Sierra Nevada, Spain, where a simple network with sub-parallel trunk streams flows down the main topographic gradient. The median lengthwise asymmetry, Δℒ, over all paired flow pathways in this mountain range is 0.19. The map projection is UTM zone 30S. Panel (b) displays an example from Sierra Madre del Sur, Mexico, showcasing a complex drainage network with tortuous flow paths. The median Δℒ over all paired flow paths in this mountain range is 1.00. The map projection is UTM zone 14Q. (c) A schematic basin with the same components as in panels (a) and (b). L_ij and L_ji, used for calculating the Δℒ of paired flow paths, are measured along the two paired flow paths (light blue) from channel heads ch_i and ch_j (yellow squares), respectively, to their common junction (red square). Δχ is determined by measuring the χ values of the channel heads starting at the junction where they merge.

Differences in network planform complexity directly control the landscape's 3D topography. For the same total relief, longer and more tortuous flow paths have diverse slope aspects and shallower channel slopes, resulting in lower local reliefs (DiBiase et al., 2010) and longer channel segments per elevation range and associated ecoclimatic zone within individual main basins. Conversely, shorter and simpler flow paths that conform to the regional gradient feature a narrower distribution of slope aspects and greater local fluvial relief, such that each main basin is expected to have shorter channel segments within any given elevation range. These characteristics affect water runoff, sediment transport capacity, rate and pattern of erosion, and the distribution of ecological niches (Rodríguez-Iturbe and Valdés, 1979; Whipple and Tucker, 2002; Badgley et al., 2017; Pelletier et al., 2018; Khosh Bin Ghomash et al., 2019; Beeson et al., 2021; Stokes and Perron, 2020).

Some of the variability in network complexity could be attributed to the level of heterogeneity in the environmental and boundary conditions affecting the landscape. Spatial gradients in tectonics (Castelltort et al., 2012; Goren et al., 2015, 2014; Habousha et al., 2023; Cowie et al., 2006; Braun et al., 2013; Mudd et al., 2022), climate (Caylor et al., 2005; Thomas et al., 2011; Abed-Elmdoust et al., 2016), and lithology (including fabric and fracture density) (Strong et al., 2019; Ward, 2019; Mudd et al., 2022) and discrete geologic structures (Hamawi et al., 2022; Scott and Wohl, 2019) are likely linked to more complex network geometry (e.g., Abed-Elmdoust et al., 2016). However, numerical studies of landscape evolution (Shelef and Hilley, 2014; Tucker and Whipple, 2002; Howard, 1994; Rinaldo et al., 1992; Sun et al., 1994b; Howard, 1990) show that variabilities in complexity emerge even when environmental and boundary conditions are spatially uniform. This means that drainage complexity could emerge from autogenic network dynamics and be independent of any heterogeneity in the applied forcings.

The same modeling studies (Shelef and Hilley, 2014; Tucker and Whipple, 2002; Howard, 1994; Sun et al., 1994b; Howard, 1990) found that changing the channel concavity index leads to variations in numerical network complexity, where drainages that are characterized by a higher concavity index are more complex. However, a similar relation was not reported in natural drainage networks, and the reasoning behind it remained elusive. The concavity index, θ, emerges as the exponent of the globally documented empirical power law relation between the drainage area, A, and slope, S, known as Flint's law (Flint, 1974; Howard, 1971; Whipple and Tucker, 1999; Willgoose et al., 1991):

$$\begin{array}{}\text{(1)}& S=K_{\mathrm{s}}{A}^{-\mathit{\theta }},\end{array}$$

where K_s is referred to as the steepness index. The concavity index describes changes in the slope of the river channel along its longitudinal profile as it accumulates drainage area downstream. Empirical studies have found that the concavity index ranges between 0.1–1, with values between 0.3–0.7 being more common (Tucker and Whipple, 2002). A sub-linear profile is characterized by concavity values close to zero, where the channel slope is nearly independent of the drainage area. In contrast, high θ values, closer to 1, indicate that most of the elevation gain occurs at higher elevations and small drainage areas (Whipple and Tucker, 1999). Channel concavity has been shown to vary with formative processes, which are primarily influenced by hydrologic conditions, i.e., the relationship between precipitation, discharge, and channel width (Whipple and Tucker, 1999; Stock and Dietrich, 2006), and by spatial gradients in tectonic uplift (Seybold et al., 2021) and climatic conditions (Roe et al., 2002). Notably, several recent studies identified links between channel concavity and prevailing climatic conditions, particularly aridity. These studies have consistently indicated that arid regions tend to exhibit lower concavity indices (Zaprowski et al., 2005; Chen et al., 2019; Getraer and Maloof, 2021; Michaelides et al., 2022).

The relationship between channel concavity and network complexity is intriguing because it suggests that the concavity index, θ, which characterizes the channel longitudinal profiles, controls the planform properties of entire branching drainage networks. One potential way in which θ may affect the planform geometry of drainage basins is through its hypothesized influence on junction branching angles (Howard, 1971, 1990; Sólyom and Tucker, 2007; Hooshyar et al., 2017; Strong and Mudd, 2022). Larger θ values are associated with larger branching angles, potentially leading to wider basins, while lower θ values are associated with smaller junction angles and narrower basins. However, consistent changes in the local metric of junction branching angle, despite their potential effect on basin scaling (Yi et al., 2018), do not necessarily correspond to variations in drainage complexity.

To proceed, it is essential to establish a formal definition of drainage network complexity. The term complexity has been used in association with drainage networks referring to various metrics and geometric properties, including, for example, channel branching angle (Devauchelle et al., 2012) and Horton–Strahler order (De Bartolo et al., 2016). A formal mathematical definition was proposed by Ranjbar et al. (2020), who used an entropy measure applied to series that describe the width and incremental area functions within a basin (Ranjbar et al., 2020, and references therein). These functions capture the variability in drainage density (channel pixels) and incremental contributing area along the catchment, from the furthest drainage divide to the outlet. However, by reducing the two-dimensional geometry of the drainage network to one-dimensional functions (Gangodagamage et al., 2014), these measures cannot account for internal basin asymmetry and the overall degree of channel tortuosity.

Ranjbar et al. (2020) further identified a significant correlation between their complexity measure and the drainage network topology, defined by the c-parameter of the self-similar Tokunaga tree constructed based on the drainage network (Pelletier and Turcotte, 2000; Zanardo et al., 2013) and describing the properties of side branching. This suggests that a topological description based on the Tokunaga tree might also be informative for network complexity. Nevertheless, the Tokunaga tree properties rely on average side branching by order, which overlooks internal basin asymmetry. Additionally, the Tokunaga-derived c-parameter assumes topologically self-similar basins, a potentially restrictive assumption.

Here, we adopt an alternative approach for defining and quantifying network complexity based on the lengthwise asymmetry between paired flow pathways that diverge from a single divide and rejoin at a junction or a common base level following Shelef and Hilley (2014) (Fig. 1c). The link between lengthwise asymmetry, capturing the multi-scale network geometry, and drainage complexity can be understood by considering two end-members. Extremely tortuous and complex networks are linked to large lengthwise asymmetry (Fig. 1b), as one tortuous flow pathway can be meaningfully longer than its across-divide pair. Conversely, simple geometry is linked to small asymmetry (Fig. 1a). Following this definition, we explore the associations between the concavity index and the network complexity as reflected by lengthwise asymmetry. We target observations from natural drainage networks and process-based rationale to quantify and better understand a climate-dependent, first-order control on the 3D geometry of landscapes and the river networks that drain them.

2 Quantifying the complexity and stability of branching drainage networks

Drainage divides delineate drainage basins, and, consequently, the planform geometry of a drainage network is closely tied to the associated drainage divide network (Shelef, 2018; Scherler and Schwanghart, 2020b; Habousha et al., 2023). The stability of the drainage network planform geometry is, therefore, tied to the stability of the divide network. When the divides shift or jump, the drainage basin's geometry changes and the network geometry is not steady. In contrast, as long as the divides remain stationary, the planform geometry of the drainage system is stable and the network topology is fixed. The stability of drainage divides can be assessed using the gradient of the parameter χ in between channel heads across divides (Willett et al., 2014), where χ is proportional to the expected steady-state elevation and is defined as (Perron and Royden, 2013)

$$\begin{array}{}\text{(2)}& \mathit{\chi }\left(x\right)=\underset{x_{\mathrm{b}}}{\overset{x}{\int }}{\left({\displaystyle \frac{A_{\mathrm{0}}}{A\left(x^{\prime }\right)}}\right)}^{\mathit{\theta }}\mathrm{d}{x}^{\prime }.\end{array}$$

In Eq. (2), x is the spatial coordinates measured upstream the channel, x_b represents the base level, x^′ is the integration parameter, and A₀ is a reference drainage area introduced to ensure that the dimension of χ [L] is independent of the value of θ.

Considering the χ values of channel heads (Fig. 1c), when the environmental conditions are spatially uniform, a zero or sufficiently small (Shelef and Goren, 2021) Δχ across a divide indicates that the divide is stable (Willett et al., 2014; Shelef and Hilley, 2014), while a large Δχ across a divide could indicate a migrating divide (Willett et al., 2014; Beeson et al., 2017; Habousha et al., 2023). To compare the stability of divides and drainages of different scales, a normalized χ difference (referred to as χ difference) across a divide is defined for paired flow pathways, which originate from two channel heads, i and j, across a single divide that join at a downstream junction or base level (Fig. 1c):

$$\begin{array}{}\text{(3)}& \mathrm{\Delta }{\mathit{\chi }}_{ij}=\mathrm{\Delta }{\mathit{\chi }}_{ji}={\displaystyle \frac{\mathrm{2}|{\mathit{\chi }}_{ij}-{\mathit{\chi }}_{ji}|}{{\mathit{\chi }}_{ij}+{\mathit{\chi }}_{ji}}},\end{array}$$

where χ_ij (χ_ji) is the χ value at channel head i (j), where χ is computed between the channel head and the common junction of flow pathways i and j.

A formal quantification of network complexity is defined in a similar manner as a measure of normalized lengthwise difference (referred to also as asymmetry) between paired flow pathways (Fig. 1c):

$$\begin{array}{}\text{(4)}& \mathrm{\Delta }{\mathcal{L}}_{ij}=\mathrm{\Delta }{\mathcal{L}}_{ji}={\displaystyle \frac{\mathrm{2}|L_{ij}-L_{ji}|}{L_{ij}+L_{ji}}},\end{array}$$

where L_ij and L_ji are the along-flow distances from the two channel heads to their common junction or base level.

These definitions offer valuable insights into how θ might influence network complexity, Δℒ. In the extreme case of θ=0, χ_ij=L_ij, and χ differences across the divide reduces to lengthwise asymmetry (Shelef and Hilley, 2014). In this case, stable planform configurations with Δχ_ij=0 for all i and j reduce to Δℒ_ij=0 for all i and j. Consequently, for θ=0, the only stable branching network is one with perfect lengthwise symmetry, the simplest possible network. As θ increases, the drainage area distribution along the flow pathways, Eq. (2), plays a growing role, such that equal χ across divides could also be achieved when Δℒ_ij≠0, as long as the drainage area distribution compensates for the lengthwise asymmetry.

3 Methods

We explore correlations between landscapes' channel concavity indices and their fluvial branching network complexity while accounting for the networks' stability. The analysis targets natural drainage networks, numerical networks generated using the surface process model DAC (Goren et al., 2014), and numerical optimal channel network simulations (Rinaldo et al., 1992).

3.1 Elongated natural mountain ranges

To explore the correlation between channel concavity and drainage network complexity in natural fluvial drainage networks, we independently quantify θ, Δℒ, and Δχ along 18 elongated mountain ranges across the globe. We choose to focus on elongated ranges (rather than study general networks) because (i) such ranges represent topographic units, whose base-level boundaries are relatively well defined; (ii) each of the ranges is relatively simple in terms of its tectonic setting, where the main faults bound the range rather than transect it; and (iii) for each range, θ and Δℒ are quantified over a relatively large domain with an along-range length between 10s–100s km. We further note that the ranges we choose are situated in both extensional and compressional settings, and, accordingly, some are bounded by normal faults and others are bounded by reverse faults. Detailed information about the elongated ranges is listed in Appendix A.

The selection of elongated mountain ranges for the current analysis adhered to specific criteria: (i) the elongated range has a single main divide from which basins drain to two opposite base levels. (ii) There should be a minimum of four basins on each flank, with the basins' outlets determined by a common elevation contour surrounding the range. (iii) The range should be free from prevalent volcanic characteristics or systematic structural control on the internal drainage pattern.

The analysis of the natural elongated mountain ranges is based on the SRTM 3 arcsec digital elevation model (DEM) (NASA Shuttle Radar Topography Mission, 2013) using the TopoToolbox topography analysis package (Schwanghart and Scherler, 2014; Scherler and Schwanghart, 2020a). The boundaries of each range were defined based on a minimum elevation threshold, chosen visually to eliminate alluvial fans and focus on bedrock rivers. Then, interstitial basins were excluded from the analyzed area, such that the analysis was based only on the main basins, draining the main divide to the boundary contour. The drainage network was extracted based on a predefined drainage area threshold, and we explore the sensitivity of the results to the drainage area threshold in Appendix B.

A single, best-suited concavity index, θ, was determined for each elongated mountain range based on the extracted drainage network and using the disorder scheme (Gailleton et al., 2021; Hergarten et al., 2016; Mudd et al., 2018). The disorder scheme assesses the extent to which the elevation-based order of channel pixels aligns with their χ-value order. The scheme involves calculating a normalized measure of pixel disorder, D^∗(θ), for predefined values of θ (Gailleton et al., 2021). The most probable θ value (referred to as the best-suited θ) is the one that minimizes the D^∗(θ) metric. For each value of θ, the calculation involves calculating $R\left(\mathit{\theta }\right)={\sum }_{i=\mathrm{1}}^n\left|{\mathit{\chi }}_{i+\mathrm{1}}\right(\mathit{\theta })-{\mathit{\chi }}_i(\mathit{\theta }\left)\right|$ when χ is sorted by elevation. Then, because the absolute value of χ depends on θ, R(θ) is normalized by χ_max(θ) to define $D\left(\mathit{\theta }\right)=\left(R\right(\mathit{\theta })-{\mathit{\chi }}_{\text{max}}(\mathit{\theta }\left)\right)/{\mathit{\chi }}_{\text{max}}\left(\mathit{\theta }\right)$ (Hergarten et al., 2016; Gailleton et al., 2021). Finally, D^∗(θ) is defined as $D^{\ast }\left(\mathit{\theta }\right)=D\left(\mathit{\theta }\right)/D_{\text{max}}$ (Gailleton et al., 2021).

The uncertainty in θ is derived based on the uncertainty in D^∗ following Gailleton et al. (2021). To evaluate it, a set of D^∗(θ) values was generated through bootstrapping iterations, with the number of iterations being 1.5 times the number of main basins in each range. In each iteration, D^∗(θ) was computed based on the drainage network of a random selection of 90 % of the main basins. The specific parameters used in the bootstrapping iterations were chosen heuristically for their relatively consistent results. The identification of the best-suited θ using the disorder scheme does not assume any specific functional dependency between the elevations and χ values but does assume that the parameters affecting the elevation along the drainage network, such as tectonic, climate, and lithology, might vary as a function of χ.

Δℒ and Δχ were computed for all divide points with distance from divide endpoints exceeding 1000 m (Scherler and Schwanghart, 2020a). For each of these divide points, Δℒ and Δχ were calculated based on the L and χ values of the two opposing (across the divide point) nearest drainage network pixels along the D8 flow routing raster and the L and χ values of the junction (or base level) of the two opposing pathways. To eliminate the influence of overall range asymmetry (which could stem from orographic effects on climate, tectonic advection, or tectonic tilting), the main divide points from which flow diverges to the two opposite base levels of the elongated range were excluded from the analysis. The calculation of χ values employed the best-suited θ value of the range. In cases where the two nearest network pixels were not at the same elevation, a correction was applied to Δℒ and Δχ. This correction adds the values of $A_{\mathrm{0}}^{\mathit{\theta }}\mathrm{\Delta }{z}_{ij}/K_{\mathrm{s}}$ and $\mathrm{\Delta }{z}_{ij}/K_{\mathrm{s}}{A}_i^{-\mathit{\theta }}$ to the χ and L of the lower pixel, respectively, thus estimating the χ and L values for channel heads that are at exactly the same elevation. Here, Δz_ij represents the elevation difference between the two nearest network pixels, K_s is the best-fit steepness index (derived from the slope of the χ–z data of the range), and A_i denotes the drainage area at the lower nearest network pixel (assuming it is labeled as i).

3.2 DAC simulations

The DAC landscape evolution model is a process-based model presented in Goren et al. (2014). DAC implements an implicit solver of the stream power incision model (Howard, 1994; Whipple and Tucker, 1999), E=K(PA)^mSⁿ, where E [L/T] is the erosion rate, K [L^(1−3m)/T^(1−m)] is the erodibility coefficient, P [L/T] is the precipitation rate, A [L²] is the drainage area, S [L/L] is the channel gradient, and m and n are positive exponents. Upon identifying $K_{\mathrm{s}}=(E/K{P}^m{)}^{\mathrm{1}/n}$ and $\mathit{\theta }=m/n$, the stream power model can be shown to reduce to Eq. (1). The solver is built upon a triangular, sparse, dynamically adjusting grid. Unlike previous implementations (Goren et al., 2014, 2015; Habousha et al., 2023) that solved for the divide location and identified captures following divide breaching, the DAC implementation used here assumes a strict steepest descent algorithm for flow routing. This choice allows better preservation of the initial conditions (starting with a random subdued topography between 0–1 m) and facilitates the comparison of drainage network complexity across values of the concavity index.

For the current analysis, we ran simulations over a domain size of 200 km × 60 km, producing elongated numerical mountain ranges with a single main water divide. The simulations apply a precipitation rate of P = 1 m yr⁻¹; an uplift rate of U = 0.5 mm yr⁻¹; and a slope exponent, n, of 1. The area exponent, m, was varied between the simulations and was maintained spatially uniformly in each simulation. The applied concavity index was calculated as $\mathit{\theta }=m/n$. The drainage density, a function of grid spacing, is independent of the model parameters and similar across the simulations. To maintain a consistent global relief despite the changes in θ, the erodibility coefficient K is adjusted across the simulations. Table C1 lists the values of K.

The simulations were run for 100 million years, ensuring topological stability by verifying that no alterations in flow routing occurred during the final 10 million years of each simulation. To ensure that the observed effects are indeed related to θ and not influenced by K, Fig. C1 replicates a subset of the analysis while keeping K constant.

In the DAC simulations, the calculation Δℒ and Δχ across divides accounts for all grid-based channel head pairs that share a divide. When a pair of channel heads does not have the same elevation, a similar correction to the one described for the natural elongated mountain ranges was used with the applied steepness index, K_s.

3.3 Optimal channel network simulations

Optimal channel network (OCN) theory (Rodriguez-Iturbe and Rinaldo, 2001; Rinaldo et al., 1992; Molnár and Ramírez, 1998; Banavar et al., 2001) suggests that natural drainage networks self-organize in a way that minimizes global energy expenditure during water flow down the network. In this general view of landscape organization, the energy is evaluated based on the network's geometry and topology, represented as sets of nodes and edges, where edges represent the channel connections between nodes. Each internal node has a unique path leading to an outlet node, and the flow paths are loopless. The total energy expenditure, denoted as P, for any network defined over the node set is determined by the sum of the local energy expenditures, P_i, along each edge, i (Sun et al., 1994b):

$$\begin{array}{}\text{(5)}& P=\sum _i{P}_i\propto \sum _i{Q}_i{S}_i{l}_i\propto \sum _i{A}_i^{\mathit{\gamma }}{l}_i=P_{\text{eq}}.\end{array}$$

S_i is the slope across edge i, A_i is the upstream drainage area (proxy for discharge Q) of node i from which edge i originates, and l_i is the length of edge i. The term on the right-hand side of Eq. (5) is referred to as the energy equivalent, P_eq. The area exponent, γ, is expected to correlate inversely with θ, the concavity index (Eq. 1) (Strong and Mudd, 2022). However, the interdependence of these two exponents as a function of environmental conditions and network hydrology, and their consequent functional relation, remains debated (Strong and Mudd, 2022). Here, we follow the formulation of Sun et al. (1994b) and define $\mathit{\gamma }=\mathrm{1}-\mathit{\theta }$.

Natural drainage networks were found to resemble numerically generated networks in a state of a local energy minimum (Rodriguez-Iturbe and Rinaldo, 2001; Colaiori et al., 1997). A commonly used criterion to identify and construct such networks is to ensure that the total energy (i.e., Eq. 5) is not reduced by a single-edge flip. An edge flip is an operation that redirects an edge that emerges from node i and used to end at node j, an immediate neighbor of i, to one of its other immediate neighbors, k≠j, without creating loops. A network configuration where any edge flip will only increase its total energy content defines a local energy minimum and corresponds to a topologically stable configuration.

We perform simulations using an iterative greedy algorithm following Rodríguez-Iturbe et al. (1992) to explore the effect of θ on networks' evolution toward the local energy minimum and the complexity of the emerging stable networks. The algorithm starts from a random network configuration, attempts a random edge flip in each iteration, and accepts the new configuration only if the edge flip reduces the total energy. Note that this approach differs from the simulated annealing algorithm (e.g., Sun et al., 1994a, b) that defines a temperature-dependent probability to accept edge flips that increase the total energy as a means to exit local minima and identify a global minimum configuration. Here, we use the greedy approach to eliminate probability-dependent changes that obscure the topological relation between stable networks generated with different θ values from the same initial random state (Rodríguez-Iturbe et al., 1992).

Optimizations with different θ are initiated from the same random network with a domain size of 200 over 60 nodes. Each node in the domain can drain to one of its eight neighbors, leading edges to be longer in the diagonal than in the rook directions. The nodes along the domain boundary are defined as outlets. Each simulation performed an optimization with a predefined θ value. The OCN approach lacks a hillslope domain, resulting in a high drainage density. Therefore, to avoid signal overwhelming by channel head pairs that merge at a single node downstream, Δℒ and Δχ are calculated only for neighboring channel head pairs that drain to different outlets (boundary nodes). Each simulation was run for 3.6 × 10⁶ iterations.

Figure 2 Complexity and concavity index. (a) Relations between lengthwise asymmetry, Δℒ (blue) and Δχ (green), and the concavity index θ. The data represent 18 elongated mountain ranges, visualized as circles on the right-hand side of the map. The circles are color-coded by the logarithm of the aridity index (AI) (Zomer et al., 2022), with yellower circles corresponding to more arid conditions. (b) Relations between Δℒ, Δχ, and the concavity index θ for numerical ranges from the DAC process model (Goren et al., 2014). (c) Relations between Δℒ, Δχ, and the concavity index θ for numerical ranges derived from simulations using an optimal channel network model. In all panels, the squares display median values and the vertical error bars indicate the 25th and 75th percentiles. The horizontal error bars in panel (a) represent the uncertainty in θ, where the square is located at the best-suited θ value. A relatively high regression slope and significant correlation is observed between θ and Δℒ and between θ and the spread of Δℒ in the natural and numerical ranges. A weaker (a, c) or nonexistent (b) correlation is observed between Δχ and θ.

4 Results – concavity correlates with lengthwise asymmetry in natural and numerical mountain ranges

Among the 18 natural elongated mountain ranges studied, a larger θ correlates with a higher median Δℒ (Fig. 2a; blue squares) and a wider spread of Δℒ values (Fig. 2a; blue bars, denoting the 25th and 75th percentiles of the Δℒ distribution in each range). Though the correlation is not necessarily linear, we quantify it through Pearson's linear correlation. The correlation coefficient between the best-suited θ and the median Δℒ is 0.92 (with a slope of 2.16 and a P-value of 4.22 × 10⁻⁸), and the correlation coefficient between the best-suited θ and the difference between the 75th and 25th percentiles of the Δℒ distribution is 0.81 (with a slope of 1.30 and a P-value of 5.54 × 10⁻⁵). These trends indicate that those natural networks that are characterized by a higher value of θ are more complex (with a larger median Δℒ) and show greater variability in their level of complexity. In contrast, low-θ natural networks have lower complexity and complexity variation.

To explore the relation between planform stability and θ, Fig. 2a also shows the correlation between Δχ and θ (green symbols and bars). The Pearson's linear correlation coefficient between the best-suited θ and the median Δχ is 0.58 (with a slope of 0.33 and a P-value of 0.01). The correlation coefficient relating θ to the difference between the 75th and 25th percentiles of the Δχ distribution is 0.61 (with a slope of 0.43 and a P-value of 0.007). Therefore, the association and sensitivity (i.e., correlation and slope) between Δχ and its spread to θ are weaker than those between Δℒ and its spread and θ. Consequently, the degree of network instability, as quantified by Δχ, cannot be invoked as a main driver of the variability in complexity for the analyzed ranges.

We cannot fully exclude the possibility that heterogeneity in the environmental conditions across the analyzed mountain ranges affect both θ and Δℒ. To address this possibility, we apply a similar analysis over two types of synthetic, steady-state landscapes of uniform environmental conditions generated using (i) the DAC process-based landscape evolution model (Goren et al., 2014) and (ii) the process-independent greedy OCN model. In both models, we ran simulations with pre-defined θ values and measured Δχ and Δℒ over the emerging drainage networks after they achieved topological stability. For Δℒ, models' results, Fig. 2b and c (blue), show similar trends to those documented for the natural elongated mountain ranges. The median and the spread of Δℒ increase with increasing θ. For Δχ (blue), model results show constant and infinitesimal values, independent of θ in the DAC simulations and with a weak dependency on θ in the OCN simulations. Notably, while the range of Δℒ has the same order of magnitude in the natural elongated ranges (Fig. 2a) and the numerical drainages (Fig. 2b and c), the Δℒ values in the simulations are mostly smaller than the natural Δℒ for the higher θ values, potentially revealing the effect of environmental heterogeneity on the natural terrains' complexity, consistent with the higher values of Δχ for these ranges.

Figure 3 Changes in lengthwise asymmetry, Δℒ as a function of concavity index, θ with differing Hack's law parameters; Eq. (6). Δℒ is calculated for idealized sub-basins along paired flow pathways, i and j, sharing a stable divide (Fig. 1c). (a) The color map represents the logarithm of Hack's coefficient ratio, assuming that Hack's exponents are identical, $h_i=h_j=\mathrm{2}$. (b) The color map represents the difference between Hack's exponents, assuming Hack's coefficients are identical and h_i=2. In both panels, L_i = 50 km, L_j≤L_i, and x_c = 0.25 km. The black and gray curves represent fitted power laws to the Δℒ–θ trends based on the natural elongated mountain ranges shown in Fig. 2a. For the median Δℒ, the fit is Δℒ=2.12θ^1.60 (black curves); for the 75th percentile, the fit is Δℒ=2.49θ^1.23 (upper gray curves); for the 25th percentile, the fit is Δℒ=1.54θ^2.14 (lower gray curves).

5 Discussion

5.1 Hack's law explains the relation between θ and Δℒ

To explain the observed correlation between θ and landscape complexity as quantified by Δℒ, we turn to an analysis of idealized channels. We set two coordinate systems, x, that follow paired flow paths from their common junction or base level to their common divide, such that x=0 is at the junction and x=L_i (x=L_j) is the common divide when measured along channel i (j) (Fig. 1c). We assume that the channels obey Hack's law (Rigon et al., 1996), such that the drainage area distribution along the sub-basins whose outlet is the common junction are expressed as follows:

$$\begin{array}{}\text{(6)}& \begin{array}{rl}{A}_i\left(x\right)& =k_{a_i}(L_i-x{)}^{h_i}\text{ for }\mathrm{0}\le x\le L_i-x_{\mathrm{c}}=L_{ij}\\ A_j\left(x\right)& =k_{a_j}(L_j-x{)}^{h_j}\text{ for }\mathrm{0}\le x\le L_j-x_{\mathrm{c}}=L_{ji},\end{array}\end{array}$$

where $k_{a_i}$ and $k_{a_j}$ are Hack's coefficients and h_i and h_j are Hack's exponents. The hillslope length, x_c, measured between the divide and the channel heads, is assumed to be uniform. We further assume that the channels obey a power law relation between the slope and the drainage area with a concavity index θ, Eq. (1), and, consequently, that the χ values at the channel heads could be defined by combining Eqs. (6) and (2):

$$\begin{array}{}\text{(7)}& {\displaystyle }\begin{array}{rl}& {\mathit{\chi }}_i(L_i-x_{\mathrm{c}})=\underset{x^{\prime }=\mathrm{0}}{\overset{x^{\prime }=L_i-x_{\mathrm{c}}}{\int }}{\displaystyle \frac{A_{\mathrm{0}}^{\mathit{\theta }}\mathrm{d}{x}^{\prime }}{A_i(x^{\prime }{)}^{\mathit{\theta }}}}\\ & =\left\{\begin{array}{ll}\frac{A_{\mathrm{0}}^{\mathit{\theta }}}{(\mathrm{1}-h_i\mathit{\theta })k_{a_i}^{\mathit{\theta }}}\left(L_i^{\mathrm{1}-h_i\mathit{\theta }}-x_{\mathrm{c}}^{\mathrm{1}-h_i\mathit{\theta }}\right)& \text{for }{h}_i\mathit{\theta }\ne \mathrm{1}\\ \frac{A_{\mathrm{0}}^{\mathit{\theta }}}{k_{a_i}^{\mathit{\theta }}}\ln \left(\frac{L_i}{x_{\mathrm{c}}}\right)& \text{ for }{h}_i\mathit{\theta }=\mathrm{1},\end{array}\right.\end{array}\text{(8)}& {\displaystyle }\begin{array}{rl}& {\mathit{\chi }}_j(L_j-x_{\mathrm{c}})=\underset{x^{\prime }=\mathrm{0}}{\overset{x^{\prime }=L_j-x_{\mathrm{c}}}{\int }}{\displaystyle \frac{A_{\mathrm{0}}^{\mathit{\theta }}\mathrm{d}{x}^{\prime }}{A_j(x{)}^{\mathit{\theta }}}}\\ & =\left\{\begin{array}{ll}\frac{A_{\mathrm{0}}^{\mathit{\theta }}}{(\mathrm{1}-h_j\mathit{\theta })k_{a_j}^{\mathit{\theta }}}\left(L_j^{\mathrm{1}-h_j\mathit{\theta }}-x_{\mathrm{c}}^{\mathrm{1}-h_j\mathit{\theta }}\right)& \text{ for }{h}_j\mathit{\theta }\ne \mathrm{1}\\ \frac{A_{\mathrm{0}}^{\mathit{\theta }}}{k_{a_j}^{\mathit{\theta }}}\ln \left(\frac{L_j}{x_{\mathrm{c}}}\right)& \text{ for }{h}_j\mathit{\theta }=\mathrm{1}.\end{array}\right.\end{array}\end{array}$$

Requiring the channels to be in a topological steady state (stable divide and stable planform configuration), the channel head χ values across the divide must be equal. For generality and simplicity, we consider the case where h_iθ≠1 and h_jθ≠1 and write the divide stability criterion, equating χ_i(L_i−x_c) to χ_j(L_j−x_c):

$$\begin{array}{}\text{(9)}& \begin{array}{rl}& {\displaystyle \frac{\mathrm{1}}{(\mathrm{1}-h_i\mathit{\theta })k_{a_i}^{\mathit{\theta }}}}\left(L_i^{\mathrm{1}-h_i\mathit{\theta }}-x_{\mathrm{c}}^{\mathrm{1}-h_i\mathit{\theta }}\right)\\ & ={\displaystyle \frac{\mathrm{1}}{(\mathrm{1}-h_j\mathit{\theta })k_{a_j}^{\mathit{\theta }}}}\left(L_j^{\mathrm{1}-h_j\mathit{\theta }}-x_{\mathrm{c}}^{\mathrm{1}-h_j\mathit{\theta }}\right).\end{array}\end{array}$$

If Hack's coefficients and exponents are identical for the two sub-basins, i.e., $k_{a_i}=k_{a_j}$ and h_i=h_j, then the only solution to Eq. (9) is lengthwise symmetry, L_i=L_j and Δℒ=0 for all values of θ.

To explore the possibility of stable topological configurations that respect equal χ across a divide while permitting Δℒ≠0, we relax the restricting assumption of an identical Hack's coefficient or exponent. Firstly, we consider the case where $k_{a_i}/k_{a_j}$ is not necessarily 1, whereas $h_i=h_j=h$. Figure 3a shows the value of $k_{a_i}/k_{a_j}$ needed to ensure equal χ across a divide, Eq. (9), as a function of θ and Δℒ. Without the loss of generality, we assume that L_j<L_i and consider values of $k_{a_i}/k_{a_j}<\mathrm{10}$, approximately a factor of 3 larger than the range of reported natural values (Montgomery and Dietrich, 1992; Mueller, 1972; Willemin, 2000; Dodds and Rothman, 2000; Shen et al., 2017; Sassolas-Serrayet et al., 2018). Considering a fixed Δℒ, the figure shows that, for low values of θ, there are no $k_{a_i}/k_{a_j}$ values within the range for which L_i≠L_j. As θ increases, Δℒ≠0 could be achieved with high $k_{a_i}/k_{a_j}$ values. As θ further increases, the $k_{a_i}/k_{a_j}$ values ensuring stable configurations with any particular Δℒ become smaller. This analysis reveals that, for small θ values, stable topologies can be achieved only with small Δℒ, but, as θ increases, small differences in Hack's coefficients permit stable topological configuration with large lengthwise asymmetry.

Next, we consider the case where h_i≠h_j, whereas $k_{a_i}=k_{a_j}=k_a$. Here, Eq. (9) becomes independent of k_a, and, for a fixed value of h_i, the value of h_j for different Δℒ can be solved only implicitly. Figure 3b shows the value of the difference h_i−h_j as a function of θ and Δℒ for h_i=2. The difference is used rather than the ratio, as with k_a (Fig. 3a), because the h exponents vary by a factor, whereas the k_a coefficients can vary by orders of magnitude. We only consider solutions where $\mathrm{0}<h_i-h_j<\mathrm{0.5}$, again representing a difference larger by a factor of approximately 3 with respect to the range of h exponents reported for natural terrains. Here, as well, for low values of θ, a solution exists only for very small Δℒ, and, as θ increases, smaller differences in the Hack's exponents allow a stable configuration with a large lengthwise asymmetry.

This analysis reveals that the widely documented geomorphic relationships of Hack's law, Eq. (6), and Flint's law, Eq. (1), are consistent with and can explain natural and numerical observations (Fig. 2) of the relations between θ and Δℒ. More specifically, the analysis shows that, when θ is large, stable configurations with zero Δχ across divides can be achieved even when Δℒ≫0, facilitated by small variations in Hack's exponent and coefficient. When θ is small, topological stability necessitates high lengthwise symmetry (small Δℒ).

The curves that overlie Fig. 3 show the best-fit power law relation between the best-suited θ and the median (black curve) and the 25th and 75th (gray curves) percentiles of Δℒ based on the elongated mountain ranges shown in Fig. 2a. The relation between the curves and the color map that underlies them reveals that high-θ natural networks achieve large Δℒ by exploiting small variations in Hack's exponent or coefficient. In contrast, low-θ networks require greater variability in Hack's parameters to achieve a much smaller Δℒ.

5.2 Optimal channel network (OCN) perspective

While the analysis based on Hack's law, Sect. 5.1, explains the observed correlation between concavity and complexity from a geomorphic scaling relations standpoint, the OCN framework helps conceptualize this correlation based on energy considerations. Figure 4a shows the evolution of the normalized energy equivalent, $P_{\text{eq}}/P_{{\text{eq}}_{\text{init}}}$, as a function of iteration number during the energy optimization process under different values of θ and starting from the same initial random network, with equivalent energy $P_{{\text{eq}}_{\text{init}}}$. The figure shows that, as θ decreases, the normalized energy equivalent reduction is greater. Figure 4b shows the evolution of the median Δℒ during the optimization procedure, displaying a similar trend to that of the normalized energy equivalent, with a greater reduction in Δℒ with decreasing θ. Notably, while the greedy algorithm ensures a monotonous energy reduction with an increasing number of iterations (Fig. 4a), the Δℒ trends are non-monotonous (Fig. 4b). The Δℒ values of the final energy-minimum networks are depicted in Fig. 2c. Appendix D discusses the trend of Δχ through the optimization iterations.

Figure 4 Optimal channel network (OCN) dynamics as a function of concavity index. Transient response and steady-state values from OCN simulations (see Sect. 3.3 for details) showing the effect of the concavity index, θ, on the reduction trends of (a) the normalized energy equivalent, $P_{\text{eq}}/P_{{\text{eq}}_{\text{init}}}$, Eq. (5); (b) the median Δℒ; and (c) the acceptance ratio of edge flip operations that reduce the total energy, i.e., the quotient of edge flips and the total number of iterations.

Examples of network topology optimized using the greedy algorithm with different values of θ are shown in Fig. 5. High θ values result in complex, tortuous networks that do not significantly differ from the random initial conditions, consistent with the low number of accepted edge flip operations (Fig. 4c). As θ decreases, the networks become less complex and the legacy of the initial conditions is gradually erased (Fig. 5c and d). A similar behavior was recorded in surface process model simulations (Shelef and Hilley, 2014; Kwang and Parker, 2019; Howard, 1994).

Figure 5 Topologies emerging from OCN simulations with different concavity indexes. Random initial conditions (top) and final, steady optimal channel networks following the application of the greedy algorithm with different θ values. Note the great complexity and the similarity to the initial conditions of the high-θ networks relative to the simple geometry of the low-θ networks that significantly differ from the initial conditions. Edges are plotted only for drainage areas greater than five nodes, and edge width scales with the number of draining nodes.

The reduction in the minimum normalized energy equivalent and network complexity observed in Fig. 4a and b and the gradual deviation from the random initial network with decreasing θ (seen in Fig. 5) could be rationalized analytically. In the limit of θ=0 and when the edge length is uniform, the energy equivalent, Eq. (5), $\sum _i{A}_i^{\mathrm{1}-\mathit{\theta }}{l}_i$, becomes proportional to $\sum _i{L}_i$ (here, l_i is the length of a single edge and L_i represents the distance to the outlet node) (Colaiori et al., 1997). In this case, the global minimum is attained when each node drains to an outlet along the shortest path, contributing its local area to the minimal number of nodes. This ensures that, in the limit of θ→0, the emerging topology is such that each L_i is minimal and therefore equal across all divides and Δℒ→0. Such an optimal network exhibits an exceptionally simple geometry that differs significantly from any random network, explaining the low-complexity configurations that characterize low-θ OCNs (Fig. 5), the many edge flips needed to achieve such configurations, and the associated large normalized energy equivalent reduction (Figs. 4a and 2c).

In the θ→1 limit, the energy equivalent becomes independent of the drainage area. Assuming uniform l_i, the energy equivalent is a function of the number of nodes in the domain and is independent of the specific drainage configuration. Consequently, all networks, including any random initial complex network, have the same minimal energy. This, in turn, explains the low number of edge flips when θ→1, the small reduction in normalized energy equivalent, and the similarity of the final optimal OCN to the initial random network.

For $\mathrm{0}<\mathit{\theta }<\mathrm{1}$, our analysis reveals a monotonous relation between θ to the normalized energy equivalent reduction (Fig. 4a), the final Δℒ (Fig. 2c), and the acceptance ratio of edge flips (Fig. 4c). Overall, the OCN analysis indicates that θ determines the multiplicity of stable topologies. Networks with a high θ value can achieve stability across a wide range of Δℒ values, allowing many possible stable topologies. This includes the formation of complex networks with large Δℒ values, similar to random networks. In contrast, to attain stability, networks with a low θ value are restricted to simple topologies with smaller Δℒ values.

5.3 Climate aridity controls network complexity

Large data compilations of river profiles revealed that channel concavity correlates with climatic and hydrologic factors. More specifically, mean annual rainfall and rainfall intensities correlate positively with the concavity index (Zaprowski et al., 2005). Likewise, the degree of aridity as quantified by the aridity index (AI), the quotient of precipitation and evapotranspiration potential (Zomer et al., 2022), correlates with channel concavity, such that in arid regions (with a low aridity index) rivers are less concave (Chen et al., 2019; Getraer and Maloof, 2021). Combining these established relations between climate and concavity index and the correlation identified here between the concavity index and drainage complexity implies that the climatic conditions at which drainage networks develop could be encoded in their complexity and thus in the large-scale planform geometry of landscapes. Consequently, arid climates, characterized by low channel profile concavity, likely favor the development of low-complexity networks, whereas a more humid climate, characterized by high-concavity channels, is expected to result in variable complexity, including high-complexity drainage networks.

We examine the relationship between the complexity, Δℒ, and the aridity index across the elongated natural mountain ranges. The aridity index for each elongated range is calculated based on pixel statistics of the global aridity index raster (Zomer et al., 2022). A mask based on the analyzed area in each elongated range is used to extract the relevant pixels from the AI raster. The median and 25th and 75th percentiles of the pixel values within each such mask are used in the analysis.

Figure 6 shows a positive correlation between climate aridity index (AI) and network complexity for the elongated mountain ranges. The Spearman's rank correlation coefficient for this correlation is 0.59 with a P-value of 0.01, indicating, that consistent with the expectations, higher complexity tends to be associated with a higher AI, representing more humid climates. It is worth noting the significance of this correlation, considering that the correlation between the concavity index and the aridity index is found to be insignificant (see Appendix E for details) and that our dataset comprises only 18 mountain ranges. This could be seen as another support for the strong link between the network planform complexity and the formative concavity index (Eq. 1), which is expected to strongly depend on the hydrologic conditions (i.e., rainfall–discharge, and discharge–channel width relations) (Whipple and Tucker, 1999; Freund et al., 2023) and could differ from the measured concavity (Seybold et al., 2021).

Figure 6 Complexity and aridity. The relation between drainage complexity, Δℒ, and aridity index, the quotient of precipitation and evapotranspiration potential (Zomer et al., 2022), for the 18 elongated mountain ranges analyzed in Fig. 2a. The box symbols represent the median values, and the bars show the 25th and 75th percentiles. A significant correlation (P=0.01) with a Spearman's rank correlation coefficient of 0.59 indicates that more complex networks (high Δℒ) are associated with more humid climatic conditions and their corresponding hydrology.

5.4 Concavity controls planform landscape evolution

The results so far reveal that high-concavity landscapes can achieve topological stability with variable complexity, whereas the stability of low-concavity landscapes is conditioned by low complexity. These findings are consequential for landscape evolution, which we further investigate and quantify using landscape evolution simulations in DAC.

5.4.1 Changing climate

A first simulation set is designed to examine how the drainage network adjusts to changes in the concavity index, reflecting changing climatic, hydrologic, and geomorphic conditions. Previous studies have linked the concavity index to channel-forming processes, where debris flow channels typically exhibit lower concavity compared to fluvial channels (Stock and Dietrich, 2006). Therefore, a decrease in the concavity index can represent aridification or a transition to a debris-dominated landscape, while an increase in concavity may indicate a transition to a more humid climate or a fluvial-dominated regime.

Figure 7 The effect of changing concavity on network complexity. Simulation results from the DAC process model (Goren et al., 2014) showing the relation between the complexity, Δℒ, and the concavity index, θ, during a two-stage scenario. Firstly, the concavity is gradually decreased by steps of 0.1 (left-pointing arrow), and, secondly, the concavity increases in steps of 0.1 (right-pointing arrows). The network and topography of the initial θ=0.9 conditions are depicted in the brown-framed topography; the θ=0.1 landscape, corresponding to the end of the first stage and the beginning of the second stage, is shown with a black-framed topography; and the final θ=0.9 landscape is shown with a purple-framed topography. Note (i) the hysteresis response of Δℒ, showing different trends depending on the directional change in the concavity index, and (ii) the difference in topographic complexity between the two θ=0.9 maps. For the current analysis, Δℒ is calculated only for channel heads that drain to different outlets, such that the shared junction is the base level. Measuring Δℒ over this longer length scale emphasizes the hysteresis signal.

The simulation set starts with a topologically stable landscape of high concavity (θ = 0.9) and high complexity (brown frame in Fig. 7). We then gradually decrease the concavity index by steps of 0.1. Following each step, we let the landscape re-equilibrate until no further topological changes are observed. We measure the median value of Δℒ for this equilibrated landscape and then use this landscape as the initial condition for the next step. The procedure continues until reaching a low concavity value of θ = 0.1. Subsequently, we gradually increase the concavity index by 0.1, following the same re-equilibration procedure, until returning to the initial concavity value of θ = 0.9.

Figure 7 shows that, during the decreasing concavity stage, the median Δℒ gradually decreases, consistent with the results shown in Fig. 2b. However, in the increasing concavity stage, a hysteresis response is observed. The median Δℒ in the increasing concavity stage is lower than that of the same concavity value in the decreasing concavity stage and overall shows only a slight increase compared to the median Δℒ of the landscape with θ=0.1 (black frame in Fig. 7). Consequently, when the concavity index returns to its initial value of θ = 0.9 (purple frame in Fig. 7), the median Δℒ of the landscape is smaller by a factor of 2.4 with respect to the median Δℒ of the initial conditions with the same concavity index.

The dynamics depicted in Fig. 7 suggest that, when aridification or a transition to a debris-flow-dominated regime takes place, there is a significant autogenic reorganization of the drainage network towards a lower-complexity configuration. In contrast, when transitioning to a more humid climate or a fluvially dominated regime, minimal reorganization is expected, and the resulting landscape may retain the complexity of the antecedent more arid or debris-flow-dominated state. This implies that the complexity of a given landscape, as reflected by Δℒ, is influenced by the lowermost concavity experienced by this landscape. For example, a low-complexity landscape, currently located in a humid climate (and exhibiting high concavity), might suggest a formation or modification history under drier (i.e., low-concavity) conditions. Differences in past aridity can therefore be invoked to explain the variability in Δℒ with aridity (Fig. 6) and the increased variability in Δℒ with concavity in the elongated mountain ranges (Fig. 2a).

Figure 8 Preservation of antecedent conditions as a function of concavity index in a DAC (Goren et al., 2014) process model simulation series. In these simulations, drainage networks evolve primarily through the incision of spatially uniform rocks, except for a narrow, elongated slab with a higher erodibility (100 times greater), a width of 1 km, and a depth that extends from the surface to 5 km into the crust. The slab is located 15 km from the southern base level. The south-draining basins originally incise into the slab-containing rocks, and the drainage network locally aligns with the high-erodibility slab, forming longitudinal channels regardless of the concavity index (topographies at the top of the figure). Following a 5 km exhumation, the slab is fully eroded away and rivers incise into uniformly erodible strata. Drainage network response to the removal of the slab depends on the concavity index. In simulations of low concavity (topographic maps on the left; a and b), significant drainage reorganization removes the longitudinal segments (topography at the bottom left; b). In contrast, high-concavity landscapes (topographic maps to the right; d and e) preserve the antecedent longitudinal segments after the slab is eroded (topography at the bottom right; e). The central graph (c) presents the preservation ratio by counting the number of longitudinal segments at 100 Myr (after slab removal by erosion) and dividing it by the number of longitudinal segments at 10 Myr (while the high-erodibility slab is still present) for all south-draining basins.

5.4.2 Fingerprints of antecedent lithologic conditions

To explore another scenario in which the complexity records the legacy of past conditions, we focus on the effect of lithology in a second set of simulations. Here, drainage networks evolve from a subdued, random topography, similar to the simulations depicted in Fig. 2b. However, in this case, a narrow, 1 km wide slab protrudes into the surrounding rocks. The slab extends down to a depth of 5 km into the crust and is positioned midway between the center of the domain and the southern base level of the evolving range, as shown in Fig. 8. The slab's erodibility is higher by a factor of 100 compared to that of the surrounding rocks. This high-erodibility slab can be conceptualized as a fault zone containing crushed, more erodible rocks.

During the initial period of approximately 10 Myr (with an imposed uplift rate of 5 × 10⁻⁴ m yr⁻¹), the developing channel networks incise into the layers of rock protruded by the higher-erodibility slab. Once the slab is fully removed by erosional exhumation of 5 km, the networks continue to incise into rocks with a uniform erodibility.

The simulation results reveal distinct behaviors of the drainage networks over time, depending on the concavity index and the presence of a high-erodibility slab. In the first 10 Myr, south-draining channel segments favor the high-erodibility slab, resulting in channel segments that develop on top of the slab (Duvall et al., 2020), parallel to the mountain range, with an east–west orientation. However, beyond this period, after the high-erodibility slab has been eroded, the response of the drainage networks diverges based on the concavity index.

To quantitatively assess this response, we analyze the number of slab-parallel channel longitudinal segments whose orientations are within 10° of the slab's orientation, where a segment connects two neighboring numerical nodes. Figure 8 illustrates the ratio of longitudinal segments after 100 Myr (post-slab removal by erosion) to the number of longitudinal segments at 10 Myr, prior to the complete exhumation of the slab, for south-draining basins. As the concavity index increases, a greater preservation of longitudinal segments is observed, although the slab is completely eroded. In the case of a landscape with θ=0.9 (right-hand side of Fig. 8), the drainage pattern predominantly retains the legacy of the high-erodibility slab, resulting in a preservation ratio of longitudinal segments close to 1. Conversely, a landscape with θ=0.2 (left-hand side of Fig. 8) is characterized by lower preservation of longitudinal segments after the removal of the high-erodibility slab, leading to a low preservation ratio.

The slab simulation set reveals that environmental conditions that sustained a high concavity index over long timescales resulted in drainage networks that recorded the cumulative effect of spatially varying heterogeneities (e.g., lithology) even long after these heterogeneities were removed. In contrast, environmental conditions associated with a low concavity index might eliminate traces of past spatial heterogeneities.

The two simulation sets indicate that landscapes with high concavity indexes are less sensitive to variations in exogenic forcing. As a result, high-concavity-index landscapes can preserve the original drainage patterns (Kwang and Parker, 2019), their specific complexity, and the legacy of the environmental conditions in which these landscapes initially formed. Conversely, networks with low concavity indexes undergo reorganization to attain a configuration of lower complexity in response to any environmental, climatic, and tectonic changes. Consequently, they possess a limited capacity to retain the legacy of the previous environmental conditions.

6 Conclusions

The current analysis reveals that the channel concavity index, θ, reflecting the channel longitudinal profile, sets a first-order control on the planform complexity of drainage networks as quantified by the statistics of asymmetry in the length of paired flow pathways, Δℒ. The variability in concavity indices thus explains the observed variability in complexity across the globe. Specifically, θ controls the multiplicity of stable planform configurations available to a drainage network. When θ is small, the number of stable configurations is small, and they are all characterized by high lengthwise symmetry, producing simple-looking drainage networks. When θ is large, the number of stable planform configurations increases, and they include configurations with a large degree of lengthwise asymmetry, producing complex geometry. Consequently, high-θ drainages can be found in topological stable configurations that are characterized by high lengthwise asymmetry (i.e., high complexity), whereas the stability of low-θ drainages is conditioned by a smaller lengthwise asymmetry (i.e., low complexity). These findings can be theoretically explained based on an energy minimization principle or by combining two empirical power laws that are readily documented across the globe: Hack's law, Eq. (6), and Flint's law, Eq. (1), describing the relation between drainage area and, respectively, channel length and channel slope. The multiplicity of steady configurations of high-θ landscapes further means that the planform geometry of these networks more readily preserves the legacy of former conditions (Kwang and Parker, 2019).

Drainage network complexity of elongated mountain ranges correlates with the aridity index, a measure of climate dryness. The correlation emerges despite the relatively small number of natural ranges we analyzed and is intuit through the effect of climate on channel formative concavity (Whipple and Tucker, 1999). Therefore, the geometric complexity of drainage networks over entire mountain ranges records information about prevailing climatic conditions.

Appendix A: Elongated mountain ranges

Comprehensive details of the 18 elongated mountain ranges and the data utilized for generating Figs. 2a and 6 are listed in Table A1. The ranges' relief and drainage networks are depicted in Fig. A1.

Table A1 Elongated Mountain Ranges¹.

¹ Data based on threshold drainage area of 145 pixels corresponding to ≈ 1 km². Values to the left of the parentheses are the median. Values in parentheses correspond to the 25th and 75th percentiles of the population.
² Latitude, longitude, and elevation of the highest peak in the range.
³ Elevation of the minimal closed contour defining the bounds of the range.
⁴ The total drainage area of the main basins in the range, draining the main divide to the mountain front, at elevation z_th used for the analysis. The calculation of θ, Δℒ, and Δχ is based on the drainage network defined by these basins, and the calculation of the aridity index statistics is based on a mask defined by the boundaries of these basins.
⁵ Apparent steepness index based on the linear regression slope of the χ–z domain.
⁶ The aridity index is calculated based on a mask defined by the main basins, based on Zomer et al. (2022).
⁷ Number of paired flow pathways analyzed for Δℒ and Δχ.
⁸ Due to the extensive area of these ranges, the concavity index calculation and the extraction of the channel heads used for the analysis of Δℒ and Δχ were based on a trimmed drainage network after eliminating channels of Strahler order 1.

Figure A1 Relief and drainage networks of the 18 elongated mountain ranges used in the analysis of Figs. 2a and 6 in the main text. The downstream extent of the drainage networks is plotted starting at a common elevation contour.

Appendix B: θ – Δℒ insensitivity to threshold drainage area in natural mountain ranges

Channels of small drainage areas are sometimes associated with a low concavity index, reflecting the dominance of debris flows (Stock and Dietrich, 2006). Therefore, we explore the sensitivity of the results to the threshold drainage area used for network extraction. Figure 1a shows the relation between θ and Δℒ for 18 elongated mountain ranges across the globe. The figure is based on drainage network extraction using a drainage area threshold of 145 pixels, corresponding to ≈ 1 km². Figure B1 shows that, when the drainage network is defined based on different thresholds, the distributions of θ, Δℒ, and Δχ vary, but the emerging trends between these parameters are generally independent of the threshold drainage area.

Figure B1 The effect of threshold drainage area, A₀, used for the extraction of the drainage network on the relations between Δℒ or Δχ and the concavity index, θ. The choice of A₀ affects the calculated values, but the positive relation between Δℒ and θ and the less significant correlation between Δχ and θ are reproduced in all cases. The three panels (a–c) show data for the same elongated ranges, except the Taiwanese Central Mountain Range, which is omitted from the left panel because the concavity estimation with A₀≈ 0.34 km² produced a bimodal distribution of values.

Appendix C: θ – Δℒ and the erodibility in the numerical DAC simulations

The erodibility coefficients used in the DAC simulations, as shown in Fig. 2b, were varied along with the value of θ to maintain an approximately equal range relief (Willett, 2010). These values are listed in Table C1.

Table C1 Concavity index and erodibility in the DAC simulations^∗.

^∗ Parameters used in the simulations depicted in Fig. 2b.

To ensure that the erodibility itself does not control Δℒ, we repeat a subset of the simulations, where the concavity index varies, and the erodibility is kept constant, resulting in numerical ranges of drastically different relief. Figure C1 shows the same trend as in Fig. 2b in the main text, indicating that this trend is independent of changing the erodibility index.

Figure C1 Relation between Δℒ (blue) and Δχ (green) and the concavity index, θ for a subject of θ values. Here, the erodibility, K, is maintained constant; however, the relations observed in Fig. 2b remain the same.

Appendix D: The evolution of Δχ during OCN simulations

The main text shows the evolution of the normalized energy equivalent and Δℒ as functions of iteration number during the greedy OCN optimization application. Here, the behavior of Δχ is shown throughout the iteration process to highlight the parallelism between the two definitions of stable networks: achieving minimal Δχ across divides and minimizing the normalized energy equivalent. Figure D1 shows that the energy minimization process also decreases the median Δχ. However, compared to the normalized energy equivalent and the Δℒ trends (Fig. 4a and b in the main text), the Δχ decrease is more uniform across values of θ, as also seen in Fig. 2c in the main text. Furthermore, the decrease in Δχ is non-monotonous.

OCNs can be linked to drainage network evolution over topography (Banavar et al., 2001). Banavar et al. (2001) found that any edge flip resulting in a reduction in the network's total energy indicates that the flip occurred towards a steeper immediate neighbor. While such a flip aligns with a decrease in Δχ between the pre- and post-flip channel heads, the non-monotonic trends of Δχ depicted in Fig. D1 indicate that the median Δχ of the network, encompassing the entire model domain, does not necessarily decrease following an edge flip. This is likely because a single flip changes the drainage area distribution along the two affected basins, which can also change the Δχ across other channel head pairs.

Figure D1 Evolution of median Δχ during the optimal channel network (OCN) simulations (see Sect. 3 in the main text) for different values of the concavity index θ.

Appendix E: Aridity index (AI), concavity index, and complexity

Figure E1 depicts the relationship between the concavity index, θ, and the aridity index, AI, for the 18 elongated mountain ranges we studied. Although a positive trend is apparent, the Spearman's rank correlation coefficient for this correlation is 0.37, and the correlation is statistically insignificant with a P-value of 0.13. In contrast, the data shown in Fig. 6 for the relationship between complexity, quantified with Δℒ and the aridity index, are statistically significant.

Figure E1 The relation between the concavity index, θ, and the aridity index, the quotient of precipitation and evapotranspiration potential (Zomer et al., 2022), for the 18 elongated mountain ranges analyzed in Fig. 2a. The box symbols represent the median values, and the bars show the 25th and 75th percentiles. Despite the positive trend, the Spearman's rank correlation is statistically insignificant (P=0.13).

This raises the question of why the aridity index better correlates with complexity than with the concavity index, particularly since previous studies, which used larger datasets, identified a correlation between θ and climate (Zaprowski et al., 2005; Chen et al., 2019; Getraer and Maloof, 2021). A possible interpretation is that the measurement of channel concavity in natural settings and a particular point in time might encompass several environmental factors, including variations in tectonics (i.e., Seybold et al., 2021) and rock types, introducing noise into the relationship between concavity index and aridity index. In contrast, our analyses, based on Hack's law and OCN theory, demonstrate that network complexity, Δℒ, is contingent on the channel formative concavity index. This index is expected to depend on hydrologic conditions (Whipple and Tucker, 1999) and therefore to exhibit a more robust correlation with climate (assuming the relative aridity of the analyzed sites did not meaningfully change since network formation).