# Channelization Cascade

###### Abstract

The hierarchy of channel networks in landscapes displays features that are characteristic of nonequilibrium complex systems. Here we show that a sequence of increasingly complex ridge and valley networks is produced by a system of nonlinear partial differential equations serving as a minimalist landscape evolution model. This channelization cascade is reminiscent of the subsequent instabilities and vortex cascade in fluid turbulence and is described by a dimensionless number accounting for the relative role of diffusive soil creep, runoff erosion, and tectonic uplift. The regularity of the patterns, generated by geometrically simple boundary conditions, is indicative of the tendency to evolve toward optimal configurations, with anomalies similar to dislocation defects observed in pattern-forming systems. The choice of specific geomorphic transport laws and boundary conditions is shown to strongly influence the channelization cascade, the spatial organization of the topographic patterns, and their hypsometric curves, thus revealing the nonlocal and nonlinear character of the underlying dynamics.

Keywords: Ridge and valley patterns Landscape evolution Detachment limited River branching Drainage area

## 1 Introduction

The spatial distribution of ridges and valleys is one of the most striking features of a landscape. It has long fascinated the physical, geomorphological, and hydrological communities, leading to the development of a rich body of work on the statistical, theoretical, and numerical analysis of landscape organization. Early works focused on the definition of stream ordering systems for the river basin characterization (e.g., Horton (1945); Strahler (1952); Shreve (1966)) and the coupled dynamics of water and sediment transport to identify stability conditions for incipient valley formation (e.g., Smith and Bretherton (1972); Loewenherz (1991); Izumi and Parker (1995)), followed by the statistical description of river networks, including scaling laws and fractal properties of river basins (e.g., Tarboton et al. (1988); Marani et al. (1991); Rodríguez-Iturbe and Rinaldo (2001); Dodds and Rothman (2000)), the related optimality principles (e.g., Rigon et al. (1993); Rodríguez-Iturbe and Rinaldo (2001)), and stochastic models (e.g., Banavar et al. (1997); Somfai and Sander (1997); Pastor-Satorras and Rothman (1998)). These studies have shed light on the spatial organization and governing statistical laws of developed river networks and explored the linkages to other branch-forming systems (Kramer and Marder, 1992; Arneodo et al., 1992; Somfai and Sander, 1997), but have not tackled the physical origin of the underlying instabilities and feedback mechanisms acting over time in the formation of the observed ridge and valley patterns (Fowler, 2011). To this purpose, landscape evolution models have been employed for the analysis of branching river networks (Perron et al., 2008, 2012) in relation to the main erosional mechanisms acting on the topography. These works represented an important step forward in the study of spatially organized patterns of ridges and valleys. However, lacking a rigorous formulation of the drainage area equation (Gallant and Hutchinson, 2011; Bonetti et al., 2018), they resort to grid-size dependent algorithms (Schoorl et al., 2000; Pelletier, 2012) coupled to ad hoc empirical adjustments (e.g., Perron et al. (2008); Pelletier (2012)), which precluded a theoretical investigation of the underlying instabilities.

In this work, we focus on landscapes characterized by runoff erosion, expressed as a function of the specific drainage area (Bonetti et al., 2018) to obtain grid-independent solutions without the introduction of additional system parameters. The resulting system of coupled, nonlinear partial differential equations (PDEs) provides a starting point for the theoretical analysis of channel-forming instabilities and landscape self-organization. The nonlocal character of the equations makes the boundary conditions extremely important. On regular domains, simulations reveal a sequence of channel instabilities reminiscent of the laminar-to-turbulent transition (Pope, 2000; Drazin and Reid, 2004; Kundu et al., 2011). The relevant dimensionless number, given by the ratio of erosion and diffusion, is thus akin to a Reynolds number. The branching sequence towards smaller and smaller valleys until soil creep becomes dominant is similar to the turbulent cascade with large scale vortices leading to smaller ones until viscous dissipation. The formation of regular pre-fractal networks of ridges and valleys brought about by the regular boundary conditions also reveals the tendency of the system to develop towards optimal configurations suggestive of maximization principles (e.g., Rigon et al. (1993)) typical of nonequiblibrium thermodynamics (Rinaldo et al., 1996; Ozawa et al., 2003; Martyushev and Seleznev, 2006; Lorente and Bejan, 2002; Bejan and Lorente, 2011; Bejan, 2016) and complex branching systems (Arneodo et al., 1992; Bensimon et al., 1986; Sander and Somfai, 2005). Our analysis is different from recent interesting contributions on groundwater-dominated landscapes (Devauchelle et al., 2012; Yi et al., 2017), where branching and valley evolution is initiated at seepage points in the landscape.

## 2 Landscape evolution in detachment limited conditions

The time evolution of the surface elevation is described by the sediment continuity equation (Dietrich et al., 2003; Perron et al., 2008; Smith, 2010; Fowler, 2011)

(1) |

where is time, is an uplift rate, and is the total sediment flux, given by the sum of fluxes related to runoff erosion/channelized flow () and soil creep processes (). The soil creep flux is assumed to be proportional to the topographic gradient (Culling, 1960, 1963), so that , being a diffusion coefficient (here assumed to be constant in space and time). In detachment-limited (DL) conditions (Howard, 1994; Izumi and Parker, 1995; Perron et al., 2008) it is assumed that all eroded material is transported outside the model domain, so that no sediment redeposition occurs. Under these conditions, the runoff erosion term can be approximated as , where is a coefficient, is the discharge per unit width of contour line, and and are model parameters. As a result, Eq. (1) becomes

(2) |

Thus the soil creep flux results in a diffusion term which tends to smoothen the surface, while the runoff erosion component is a sink term which excavates the topography as a function of local slope and specific water flux.

The surface water flux can be modeled using the continuity equation

(3) |

where is the water height, the direction of the flow, and the rainfall rate effectively contributing to runoff production. Assuming representative steady state conditions, with constant and unitary rainfall rate and water flowing at unitary speed in the direction opposite to the landscape gradient (i.e., ), Eq. (3) reduces to an equation for the specific catchment area (Bonetti et al., 2018),

(4) |

since , , and coincide numerically. The resulting landscape evolution equation is

(5) |

and is coupled to Eq. (4). It is important to observe that has units of length and is the specific drainage area (, defined per unit width of contour line , see Bonetti et al. (2018)), while most landscape evolution models use the total drainage area (see, e.g., Rodríguez-Iturbe and Rinaldo (2001); Perron et al. (2008); Pelletier (2012); Chen et al. (2014)). The value of is generally evaluated using flow-routing algorithms (e.g., D8, D); however, the direct use of in Eq. (4) produces grid-dependent results. To correct for this, the drainage area is then often modified to account for the channel width (see, e.g., Perron et al. (2008); Pelletier (2012)), but this results in ad hoc approximations with the introduction of additional empirical parameters. Conversely, the use of avoids grid-dependence of the resulting topography (see Supporting Information, SI, for results obtained for different grid resolution values).

### 2.1 Nondimensional analysis

Eq. (5) can be non-dimensionalized to quantify the relative impact of soil creep, runoff erosion, and uplift on the landscape morphology. Using a typical length scale of the domain, , and the parameters of Eqs. (4) and (5), the following dimensionless quantities can be introduced: , , , , and . With these quantities, Eq. (5) becomes

(6) |

where

(7) |

As we will see later, this describes the tendency to form channels in a way which is closely reminiscent of the global Reynolds number in fluid mechanics. A similar quantity based on an internal length scale was used in Perron et al. (2008); this is similar to a Reynolds number based on a Taylor microscale (Tennekes and Lumley, 1972), although it was called Péclet number (Perron et al., 2008) for the purely formal resemblance of Eq. (5) to an advection diffusion equation when the actual coupling with the drainage-area equation (Eq. (4)) is overlooked.

The definition of as a function of external variables also allows us to directly explore the role of the uplift rate on the formation of ridge and valley patterns. In particular, when the slope exponent is equal to 1, the relative proportion of runoff erosion and soil creep can be seen to be independent of the uplift rate. However, if the uplift acts to increase the runoff erosion component, while for it enhances the diffusion component of the system. This results in different drainage network patterns as a function of uplift rates, as will be investigated in the following sections.

### 2.2 Simulation setup and modeling scenarios

Numerical simulations were performed on a regular square grid of size , using forward differencing in time and centered difference approximations for the spatial derivatives. The total drainage area was computed at each grid point with the algorithm using the Matlab TopoToolbox. The specific catchment area was then approximated as (Tarboton, 1997), being the grid size.

A first set of simulations was run assuming and (see model parameters in the SI, Table S1). Various simulations were performed to identify the critical values of at which the system transitions to different valley/ridge patterns. We run a number of simulations having a domain size of m and covering all the different valley/ridge patterns (i.e., regimes). Additionally, simulations were run considering different domain sizes (e.g., 50, 75, 152, and 200 m, see Table S1 in the SI) to show that the final ridge/valley pattern is embedded solely in the value of .

To assess the sensitivity of the network pattern to the slope and water flux exponents, simulations were performed changing the values of and . In particular, four sets of simulations were performed (i.e., increasing/decreasing both and ), and assuming fixed domain boundaries (i.e., down-cutting rate equal to zero). Parameters for these simulation sets are provided in the SI (Table S2).

Rectangular domains were also considered to evaluate the effect of boundary conditions (i.e., domain shape and size) on the resulting ridge and valley networks. Simulations were run assuming , (same as the first set of simulations), and over domains of fixed length in the direction (i.e., m), and varying length in the direction. Simulations were carried out for = 0.5, 1, 1.5, 2, 2.5, 4.8, 5, and 5.2, where is a shape factor defined as the ratio between the two horizontal length scales and , namely . A set of additional numerical experiments over a rectangular domain with = 5 were performed for different values of the dimensionless parameter (i.e., equal to 20, 40, 80, and 126).

## 3 Results

### 3.1 Organized ridge and valley patterns

Simulation results for scenarios with and are shown in Fig. 1. As can be seen from Eq. (7), for the dimensionless parameter is independent of the uplift rate, so that the spatial patterns of Fig. 1 are only a function of the relative proportions of the soil creep and runoff erosion components (see also the SI). For low values (i.e., ) no channels are formed and the topography evolves to a smooth surface dominated by diffusive soil creep (Fig. 1a). As the runoff erosion coefficient is increased the system progressively develops one, three, and five channels on each side of the square domain for , , and , respectively (Fig. 1b-d). When is increased above the central channels develop secondary branches, with the main central channel becoming of Strahler order two (Fig. 1e). As is further increased seven channels can be observed originating on each side of the domain, and the main central channel further branches (Fig. 1f-i). Simulations were run here up to , with the central channel developing several secondary branches and becoming of Strahler order five. Thus, for and and within the range of values explored here, nine regimes can be identified based on the number of channels forming on each side of the square domain and the maximum Strahler order (Fig. 1j).

As the resulting landscape transitions from a smooth topography to a progressively more dissected one, the shape of the hypsometric curve varies from concave up to one with a concave down portion for low elevations (Fig. 1k). In particular, channel formation (with no secondary branching) causes the hypsometric curve to progressively lower as a result of the lower altitudes observed in the topography, while maintaining a concave up profile. As secondary branches develop, the hypsometric curve transitions to a concave/convex one, with the convex portion at lower altitudes becoming more evident as is increased.

The striking regularity of the drainage and ridge patterns induced by the simple geometry of the domain is reminiscent of regular pre-fractal structures (e.g., Peano basin (Mandelbrot, 1982; Marani et al., 1991; Rodriguez-Iturbe et al., 1992; Flammini and Colaiori, 1996; Rodríguez-Iturbe and Rinaldo, 2001)) and is indicative of the fundamental role of boundary conditions for the emergence of the observed regular patterns, due to the highly non-local control introduced by the drainage area term.

### 3.2 Effect of runoff erosion laws

To investigate the effect of using different geomorphic transport laws for the runoff erosion term (i.e., changing the slope and water flux exponents), simulations were run modifying and by and , respectively. Results are shown in Fig. 2 and in the SI. When the value of is different from unity, the resulting ridge/valley patterns depend on the uplift rate , as per Eq. (7). When is increased the system displays channelization and secondary branching for higher values of (i.e., points are shifted to the right in Fig. 2a-b), with a more dissected planar geometry characterized by narrower valleys and smaller junction angles (Fig. 2c-f). A decrease in the value of makes the system develop smoother geometries, with wider valleys and the first secondary branching developing when only three channels per each side of the domain are present (see Fig. 2g-j). This results in a hypsometric curve with a more pronounced basal convexity as is increased above unity, as a consequence of the progressively more dissected topography (see SI, Fig. S3).

Conversely, a variation in results in an opposite behavior of the system. In fact, when is increased (Fig. 2k-n) the system develops secondary branching when only three channels are present on each side of the domain, with the formation of less numerous but wider valleys with higher junction angles, and a reduced basal convexity in the hypsometric curve (see SI). On the other end, a decrease in the value of the exponent results in a more dissected landscape, with narrower valleys and the development of up to nine channels per domain side before the first secondary branching occurs (Fig. 2o-r). This results in a more pronounced transition of the hypsometric curve to a concave down shape for low altitudes (see SI, Fig. S3).

### 3.3 Large rectangular domains

To assess boundary-condition effects on branching patterns we also considered large rectangular domains (so that is constructed using the distance between the longest sides). In our analogy with turbulent flows, this case corresponds to the flow of viscous fluids between parallel plates (Drazin and Reid, 2004; Kundu et al., 2011). Fig. 3 compares the drainage patterns obtained as a function of for rectangular domains of size 100 m by 500 m. Analogously to results for the square domain, for small values the soil creep component dominates resulting in an unchannelized smooth topography (Fig. 3a). Progressively increasing the runoff erosion component creates evenly spaced valleys which become narrower as the runoff erosion term increases (Fig. 3b-c), coherently with results in Perron et al. (2008). Further increasing provides progressively more dissected landscapes with the emergence of secondary branching (Fig. 3d-e). As in direct numerical simulations of turbulence large Reynolds numbers reduce the Kolmogorov microscale inducing smaller and smaller vortices, here increasing leads to finer and finer branching which becomes prohibitive from a computational standpoint.

The mean elevation profiles, computed as average elevation values along the axis and neglecting the terminal parts of the domain (i.e., white shaded areas in Fig. 3a-e) to avoid boundary effects, are shown in Fig. 3g-k. As increases (i.e., the topography becomes progressively more dissected) the mean elevation profile tends to become more uniform as a result of the more complex topography with ridges formed near the domain boundaries. Such a behavior of the mean elevation profiles for increasing is similar to the uniformization of mean velocity profiles across a pipe section as the flow transitions from laminar to turbulent with increasing Reynolds number (Kundu et al., 2011).

The effect of boundary conditions on the spatial regularity of ridge and valley patterns becomes apparent also when comparing simulations with different aspect ratios (see also the SI, Figs. S4 and S7). As can be seen in Fig. 3f, when the domain size is only modified by 20 m, while the mean elevation profile remains practically invariant (Fig. 3k-l), the resulting spatial pattern is visually different: two identical adjacent valleys are formed for (black rectangle in Fig. 3f), which are not observed for as well as (see SI, Fig. S7), suggesting that there must be some optimal domain length allowing for the formation of regular ridge and valley patterns. This is also evident from an analysis of cross-sections along the longer sides of the domain (SI, Figs. S5-S6): the lack of spatial regularity of ridges and valleys in the cross sections further suggests that a proportionality between the domain dimension and the characteristic valley spacing is needed to accommodate the formation of regular ridge/valley patterns. This results in the formation of dislocation defects, as highlighted in the example of Fig. 3f, as it is typical in nonlinear pattern-forming PDEs (Cross and Hohenberg, 1993).

## 4 Discussion and conclusions

A succession of increasingly complex networks of ridges and valleys was produced by a system of nonlinear PDEs serving as a minimalist model for landscape evolution in DL conditions. The sequence of instabilities is reminiscent of the subsequent bifurcations in fluid dynamic instabilities (Cross and Hohenberg, 1993; Drazin and Reid, 2004; Kundu et al., 2011) and can be tracked by a dimensionless number (), akin to a Reynolds number, which accounts for the relative proportions of runoff erosion, soil creep, and uplift in relation to the typical domain size. The analogy with fluid turbulence carries over to: i) the competition between runoff erosion and diffusion, which resembles the competition between viscous and inertial forces in turbulence; ii) the appearance of a minimum branching scale similar to the Kolmogorov microscale (Pope, 2000); as well as iii) the parallel between the uniformization of the mean hypsometric curves as the instabilities grow and the effect of turbulence on the mean velocity profiles. Both landscape evolution and turbulence are examples of driven nonequilibrium systems, in which a hierarchical pattern develops toward finer scales (branching vs. eddie cascade).

The channelization cascade with , the spacing between valleys, the branching angles, as well as the number of secondary branches are strongly dependent upon the choice of the geomorphic transport laws, as shown by our analysis of the nonlinearity parameters and , and reflected in the discussion in the literature about the values of such exponents (Chen et al., 2014) as well as their linkage to climate, vegetation cover, and soil properties (Montgomery et al., 2001; Lowman and Barros, 2014).

Characteristic spatial configurations were shown to emerge over both square and rectangular domains from the trade-off between diffusion and erosion. The striking regularity of the ridge and valley networks has the characteristics of regular pre-fractals (e.g., the Peano basin (Mandelbrot, 1982; Marani et al., 1991; Rodriguez-Iturbe et al., 1992; Flammini and Colaiori, 1996)), suggesting that the system likely obeys some optimization principle typical of complex systems in nonequilibrium steady states (Rigon et al., 1993; Rinaldo et al., 1996), such as the maximum entropy production hypotheses (Ozawa et al., 2003; Martyushev and Seleznev, 2006) and the constructal theory (Lorente and Bejan, 2002; Bejan and Lorente, 2011; Bejan, 2016). Particularly, the ridge and valley networks of Fig. 1 highly resemble Fig. 5 in Lorente and Bejan (2002), where optimized tree-shaped flow paths were constructed to connect one point to many points uniformly distributed over an area.

The shape of the hypsometric curve depends on the level of channelization and branching (Willgoose and Hancock, 1998) and thus on the dominant erosional mechanisms acting on the landscape (i.e., interplay between runoff erosion, soil creep, and uplift) and the various landscape properties affecting diffusion and erosion coefficients, such as soil type, vegetation cover, and climate. The results for rectangular domains highlight that shape and size of the basin also play a role in the delineation of the hypsometric curve. This is likely related to the fact that a smaller domain will have a higher proportion of its area dominated by diffusion, so that its hypsometric curve will be likely to display a less pronounced basal convexity (Willgoose and Hancock, 1998). Our focus on steady state conditions, rather than the transient dynamics, misses however the differences between the hypsometry of juvenile and old landscapes. It is likely that, during the early stages of the basin development when the drainage network is formed, the hypsometric curve will present a more pronounced basal convexity (Strahler, 1952) regardless of the value of , progressively transitioning toward its quasi-equilibrium form during the “relaxation phase” (Bonetti and Porporato, 2017). It will be interesting to explore whether such effects can be captured by quasi steady state simulations with time dependent erosion laws, as well as to investigate the similarities between such slow relaxations (e.g., Fig. 3), often towards slightly irregular configurations rather than perfectly regular networks, and the presence of defects in crystals and the amorphous configurations originating in glass transition (Debenedetti and Stillinger, 2001).

## Acknowledgments

We acknowledge support from the USDA Agricultural Research Service cooperative agreement 58-3826408-3-027, and US National Science Foundation (NSF) grants CBET-1033467, EAR-1331846, EAR-1316258, EAR-1338694, and the Duke WISeNet Grant DGE-1068871.

## References

- Arneodo et al. (1992) Arneodo, A., F. Argoul, E. Bacry, J. Muzy, and M. Tabard (1992), Golden mean arithmetic in the fractal branching of diffusion-limited aggregates, Physical review letters, 68(23), 3456.
- Banavar et al. (1997) Banavar, J. R., F. Colaiori, A. Flammini, A. Giacometti, A. Maritan, and A. Rinaldo (1997), Sculpting of a fractal river basin, Physical Review Letters, 78(23), 4522.
- Bejan (2016) Bejan, A. (2016), Advanced engineering thermodynamics, John Wiley & Sons.
- Bejan and Lorente (2011) Bejan, A., and S. Lorente (2011), The constructal law and the evolution of design in nature, Physics of life Reviews, 8(3), 209–240.
- Bensimon et al. (1986) Bensimon, D., L. P. Kadanoff, S. Liang, B. I. Shraiman, and C. Tang (1986), Viscous flows in two dimensions, Reviews of Modern Physics, 58(4), 977.
- Bonetti and Porporato (2017) Bonetti, S., and A. Porporato (2017), On the dynamic smoothing of mountains, Geophysical Research Letters, 44(11), 5531–5539.
- Bonetti et al. (2018) Bonetti, S., A. D. Bragg, and A. Porporato (2018), On the theory of drainage area for regular and non-regular points, Proc. R. Soc. A, 474(2211), 20170,693.
- Chen et al. (2014) Chen, A., J. Darbon, and J.-M. Morel (2014), Landscape evolution models: A review of their fundamental equations, Geomorphology, 219, 68–86.
- Cross and Hohenberg (1993) Cross, M. C., and P. C. Hohenberg (1993), Pattern formation outside of equilibrium, Reviews of Modern Physics, 65(3), 851.
- Culling (1960) Culling, W. E. H. (1960), Analytical theory of erosion, The Journal of Geology, 68(3), 336–344.
- Culling (1963) Culling, W. E. H. (1963), Soil creep and the development of hillside slopes, The Journal of Geology, 71(2), 127–161.
- Debenedetti and Stillinger (2001) Debenedetti, P. G., and F. H. Stillinger (2001), Supercooled liquids and the glass transition, Nature, 410(6825), 259.
- Devauchelle et al. (2012) Devauchelle, O., A. P. Petroff, H. F. Seybold, and D. H. Rothman (2012), Ramification of stream networks, Proceedings of the National Academy of Sciences, 109(51), 20,832–20,836.
- Dietrich et al. (2003) Dietrich, W. E., D. G. Bellugi, L. S. Sklar, J. D. Stock, A. M. Heimsath, and J. J. Roering (2003), Geomorphic transport laws for predicting landscape form and dynamics, Prediction in Geomorphology, pp. 103–132.
- Dodds and Rothman (2000) Dodds, P. S., and D. H. Rothman (2000), Scaling, universality, and geomorphology, Annual Review of Earth and Planetary Sciences, 28(1), 571–610.
- Drazin and Reid (2004) Drazin, P. G., and W. H. Reid (2004), Hydrodynamic Stability, Cambridge Mathematical Library, 2 ed., Cambridge University Press, doi:10.1017/CBO9780511616938.
- Flammini and Colaiori (1996) Flammini, A., and F. Colaiori (1996), Exact analysis of the Peano basin, Journal of Physics A: Mathematical and General, 29(21), 6701.
- Fowler (2011) Fowler, A. (2011), Mathematical geoscience, Springer Science & Business Media.
- Gallant and Hutchinson (2011) Gallant, J. C., and M. F. Hutchinson (2011), A differential equation for specific catchment area, Water Resources Research, 47(5).
- Horton (1945) Horton, R. E. (1945), Erosional development of streams and their drainage basins; hydrophysical approach to quantitative morphology, Geological Society of America Bulletin, 56(3), 275–370.
- Howard (1994) Howard, A. D. (1994), A detachment-limited model of drainage basin evolution, Water Resources Research, 30(7), 2261–2285.
- Izumi and Parker (1995) Izumi, N., and G. Parker (1995), Inception of channelization and drainage basin formation: upstream-driven theory, Journal of Fluid Mechanics, 283, 341–363.
- Kramer and Marder (1992) Kramer, S., and M. Marder (1992), Evolution of river networks, Phys. Rev. Lett., 68, 205–208.
- Kundu et al. (2011) Kundu, P. K., I. M. Cohen, and D. W. Dowling (2011), Fluid Mechanics 5th ed., Elsevier.
- Loewenherz (1991) Loewenherz, D. S. (1991), Stability and the initiation of channelized surface drainage: a reassessment of the short wavelength limit, Journal of Geophysical Research: Solid Earth, 96(B5), 8453–8464.
- Lorente and Bejan (2002) Lorente, S., and A. Bejan (2002), Combined flow and strength geometric optimization: internal structure in a vertical insulating wall with air cavities and prescribed strength, International Journal of Heat and Mass Transfer, 45(16), 3313–3320.
- Lowman and Barros (2014) Lowman, L. E. L., and A. P. Barros (2014), Investigating links between climate and orography in the central Andes: Coupling erosion and precipitation using a physical-statistical model, Journal of Geophysical Research: Earth Surface, 119(6), 1322–1353.
- Mandelbrot (1982) Mandelbrot, B. B. (1982), The fractal geometry of nature, vol. 1, WH freeman New York.
- Marani et al. (1991) Marani, A., R. Rigon, and A. Rinaldo (1991), A note on fractal channel networks, Water Resources Research, 27(12), 3041–3049.
- Martyushev and Seleznev (2006) Martyushev, L. M., and V. D. Seleznev (2006), Maximum entropy production principle in physics, chemistry and biology, Physics reports, 426(1), 1–45.
- Montgomery et al. (2001) Montgomery, D. R., G. Balco, and S. D. Willett (2001), Climate, tectonics, and the morphology of the Andes, Geology, 29(7), 579–582.
- Ozawa et al. (2003) Ozawa, H., A. Ohmura, R. D. Lorenz, and T. Pujol (2003), The second law of thermodynamics and the global climate system: A review of the maximum entropy production principle, Reviews of Geophysics, 41(4).
- Pastor-Satorras and Rothman (1998) Pastor-Satorras, R., and D. H. Rothman (1998), Scaling of a slope: the erosion of tilted landscapes, Journal of Statistical Physics, 93(3-4), 477–500.
- Pelletier (2012) Pelletier, J. D. (2012), Fluvial and slope-wash erosion of soil-mantled landscapes: detachment-or transport-limited?, Earth Surface Processes and Landforms, 37(1), 37–51.
- Perron et al. (2008) Perron, J. T., W. E. Dietrich, and J. W. Kirchner (2008), Control on the spacing of first-order valleys, Journal of Geophysical Research, 113, F04,016.
- Perron et al. (2012) Perron, J. T., P. W. Richardson, K. L. Ferrier, and M. Lapôtre (2012), The root of branching river networks, Nature, 492(7427), 100.
- Pope (2000) Pope, S. B. (2000), Turbulent Flows, Cambridge University Press, Cambridge, UK.
- Rigon et al. (1993) Rigon, R., A. Rinaldo, I. Rodriguez-Iturbe, R. L. Bras, and E. Ijjasz-Vasquez (1993), Optimal channel networks: a framework for the study of river basin morphology, Water Resources Research, 29(6), 1635–1646.
- Rinaldo et al. (1996) Rinaldo, A., A. Maritan, F. Colaiori, A. Flammini, R. Rigon, I. Rodriguez-Iturbe, and J. R. Banavar (1996), Thermodynamics of fractal networks, Physical Review Letters, 76(18), 3364.
- Rodríguez-Iturbe and Rinaldo (2001) Rodríguez-Iturbe, I., and A. Rinaldo (2001), Fractal river basins: chance and self-organization, Cambridge University Press.
- Rodriguez-Iturbe et al. (1992) Rodriguez-Iturbe, I., A. Rinaldo, R. Rigon, R. L. Bras, E. Ijjasz-Vasquez, and A. Marani (1992), Fractal structures as least energy patterns: The case of river networks, Geophysical Research Letters, 19(9), 889–892.
- Sander and Somfai (2005) Sander, L. M., and E. Somfai (2005), Random walks, diffusion limited aggregation in a wedge, and average conformal maps, Chaos: An Interdisciplinary Journal of Nonlinear Science, 15(2), 026,109.
- Schoorl et al. (2000) Schoorl, J. M., M. P. W. Sonneveld, and A. Veldkamp (2000), Three-dimensional landscape process modelling: the effect of DEM resolution, Earth Surface Processes and Landforms, 25(9), 1025–1034.
- Shreve (1966) Shreve, R. L. (1966), Statistical law of stream numbers, The Journal of Geology, 74(1), 17–37.
- Smith (2010) Smith, T. R. (2010), A theory for the emergence of channelized drainage, Journal of Geophysical Research: Earth Surface, 115(F2).
- Smith and Bretherton (1972) Smith, T. R., and F. P. Bretherton (1972), Stability and the conservation of mass in drainage basin evolution, Water Resources Research, 8, 1506–1529.
- Somfai and Sander (1997) Somfai, E., and L. M. Sander (1997), Scaling and river networks: A Landau theory for erosion, Physical Review E, 56(1), R5.
- Strahler (1952) Strahler, A. N. (1952), Hypsometric (area-altitude) analysis of erosional topography, Geological Society of America Bulletin, 63(11), 1117–1142.
- Tarboton (1997) Tarboton, D. G. (1997), A new method for the determination of flow directions and upslope areas in grid digital elevation models, Water Resources Research, 33(2), 309–319.
- Tarboton et al. (1988) Tarboton, D. G., R. L. Bras, and I. Rodriguez-Iturbe (1988), The fractal nature of river networks, Water Resources Research, 24(8), 1317–1322.
- Tennekes and Lumley (1972) Tennekes, H., and J. L. Lumley (1972), A first course in turbulence, MIT press.
- Willgoose and Hancock (1998) Willgoose, G., and G. Hancock (1998), Revisiting the hypsometric curve as an indicator of form and process in transport-limited catchment, Earth Surface Processes and Landforms, 23(7), 611–623.
- Yi et al. (2017) Yi, R., Y. Cohen, H. Seybold, E. Stansifer, R. McDonald, M. Mineev-Weinstein, and D. H. Rothman (2017), A free-boundary model of diffusive valley growth: theory and observation, Proc. R. Soc. A, 473(2202), 20170,159.

## Supplementary Information

### Introduction

Additional results from the landscape evolution model are presented in this Supplementary Information. In particular, a sensitivity analysis was performed to investigate the effect of grid size and downcutting rate on the drainage patterns obtained from the numerical experiments (results are shown in Figs. S1-S2 and discussed below). Hypsometric curves obtained for different values of the water flow exponent and slope exponent are shown in Fig. S3. Additional results for rectangular domains are provided in Figs. S4-S7 and discussed below. Lastly, parameter values for the first simulation set (i.e., with and ) and for simulations with varying exponents and are given in Tables S1 and S2, respectively.

### Effect of grid resolution and uplift/downcutting rate

As explained in the main text, using the specific drainage area as a proxy for the water flux ensures the modeled ridge/valley patterns to be robust to grid size variations, while theoretically justified in terms of a continuity equation. To assess the independence of such drainage patterns on the specific grid resolution chosen in the simulations, additional numerical tests were performed. In particular, simulations were run over the square domain of size m for = 40, 125, and 200 assuming grid sizes equal to 0.5 and 2 m. Results are compared to those obtained for m in Fig. S1, confirming that the resulting patterns and hypsometric curves are independent of the grid size value .

As we shall see from Eq. (7) in the main text, when the value of is independent of the uplift rate (or, analogously, the constant downcutting rate imposed at the four boundaries). To validate this statement additional simulations were run for = 40, 125, and 200 and considering downcutting rates at the boundaries of 50, 500, and 5000 m/yr. Results (Fig. S2) confirm that the geometry of ridge and valley lines as well as the normalized hypsometric curves are independent of the value of the uplift rate when .

### Results for rectangular domains

Results for rectangular domains with different shape factors are provided in Fig. S4. The drainage network shown in Fig. S4b is the result obtained over a square domain () for : in this case the system is in regime VI (see Fig. 1 in the main text) with five channels per each side of the domain and maximum Strahler order equal to two (i.e., one secondary branching). Reducing the shape factor does not provide the system with a drainage area sufficient to actually develop large valleys so that only parallel channels are formed (Fig. S4a). This is related to the fact that should be computed using as length scale the most constraining size of the domain. In fact, for this simulation the value of computed with m is 70.7, meaning that the system is in regime III (i.e., three channels per side, which is actually observed on the shorter side of the domain, Fig. S4a). Increasing above unity (Fig. S4c-f) provides patterns where the valley on the shorter side of the domain (i.e., m) still develops only one secondary branch, coherently with the expected behavior for . However, the number of channels on can be higher than five, as the valleys have more freedom to adjust with the rest of the geometry, thus partially breaking the spatial symmetry of the square case. Furthermore, on the longer side of the domain, channels tend to develop more secondary branches (i.e., Strahler order ). This is a result of the fact that, while on the square domain the two main ridgelines are necessarily formed along the diagonals thus providing a fixed and equal partition of on the four sides of the domain, the valleys developing on longer sides of rectangular domains have more space to widen and the collection of additional drainage area allows further branching.

The normalized hypsometric curves for rectangular domains of shape factor = 0.5, 2, and 5 are compared with that of the square domain (i.e., ) in Fig. S4g. For the hypsometric curve is mainly concave up (as we would expect given that there is no secondary branching), with a slight convexity for very high elevations. As is increased above unity, the progressively more dissected landscape (i.e., with multiple secondary branches) results in a more pronounced basal convexity of the normalized hypsometric curve (see, for example, results for ).

Figs. S5-S6 show cross-sections along the longer side of rectangular domains for different values of the dimensionless parameter . It is interesting to note that, while the probability distributions of the elevation profiles are similar between cross-sections on each side of the main central ridgeline (red and black, respectively), the lack of spatial synchronicity of ridges and valleys in the cross sections suggests that a proportionality between the domain dimension and the characteristic valley spacing is needed to accomodate the formation of regular ridge/valley patterns. Fig. S7a-c compares the drainage networks obtained for and shape factor equal to 4.8, 5, and 5.2 (i.e., was changed by 20 m in each numerical experiment). It is interesting to note that, while the mean elevation profiles (Fig. S7d) are practically identical, the spatial patterns that are formed are visully different. In particular, two identical adjacent valleys are formed for (black rectangle in Fig. S7b), which are not observed for and . Furthermore, dislocation defects can be observed in all simulations (blue dashed rectangles).

Regime | |||||||||

[1/(myr)] | [m/yr] | [m/yr] | [m/yr] | [m] | [m] | ||||

7.1 | 1.010 | 510 | 510 | 0.05 | 50 | 0.5 | 0 | 0 | I |

20.0 | 1.010 | 510 | 510 | 500 | 100 | 1 | 0 | 0 | I |

28.3 | 4.010 | 510 | 510 | 100 | 50 | 0.5 | 0 | 0 | I |

32.5 | 4.610 | 510 | 510 | 100 | 50 | 0.5 | 1 | 1 | II |

37.5 | 1.010 | 510 | 510 | 100 | 152 | 2 | 1 | 1 | II |

40.0 | 2.010 | 510 | 510 | 500 | 100 | 1 | 1 | 1 | II |

44.2 | 6.310 | 510 | 510 | 0.05 | 50 | 0.5 | 1 | 1 | II |

56.6 | 1.010 | 510 | 510 | 100 | 200 | 2 | 1 | 1 | II |

60.0 | 3.010 | 510 | 510 | 100 | 100 | 1 | 1 | 3 | III |

62.5 | 3.110 | 510 | 510 | 500 | 100 | 1 | 1 | 3 | III |

70.0 | 3.510 | 510 | 510 | 500 | 100 | 1 | 1 | 3 | III |

80.0 | 4.010 | 510 | 510 | 500 | 100 | 1 | 1 | 3 | III |

90.0 | 4.510 | 510 | 510 | 500 | 100 | 1 | 1 | 3 | III |

95.0 | 4.810 | 510 | 510 | 500 | 100 | 1 | 1 | 3 | III |

100.0 | 5.010 | 510 | 510 | 500 | 100 | 1 | 1 | 5 | IV |

110.0 | 5.510 | 510 | 510 | 100 | 100 | 1 | 1 | 5 | IV |

125.0 | 6.310 | 510 | 510 | 500 | 100 | 1 | 1 | 5 | IV |

141.4 | 2.010 | 510 | 510 | 0.05 | 50 | 0.5 | 1 | 5 | IV |

150.0 | 7.510 | 510 | 510 | 100 | 100 | 1 | 1 | 5 | IV |

160.0 | 8.010 | 510 | 510 | 500 | 100 | 1 | 2 | 5 | V |

180.0 | 9.010 | 510 | 510 | 500 | 100 | 1 | 2 | 5 | V |

190.0 | 9.510 | 510 | 510 | 500 | 100 | 1 | 2 | 7 | VI |

200.0 | 1.010 | 510 | 510 | 500 | 100 | 1 | 2 | 7 | VI |

220.0 | 1.110 | 510 | 510 | 500 | 100 | 1 | 2 | 7 | VI |

234.2 | 6.310 | 510 | 510 | 0.05 | 152 | 2 | 3 | 7 | VII |

240.0 | 1.210 | 510 | 510 | 500 | 100 | 1 | 3 | 7 | VII |

259.8 | 2.010 | 510 | 510 | 5 | 75 | 0.5 | 3 | 7 | VII |

272.8 | 2.110 | 510 | 510 | 500 | 75 | 0.5 | 3 | 7 | VII |

282.8 | 5.010 | 510 | 510 | 100 | 200 | 2 | 4 | 7 | VIII |

311.1 | 5.510 | 510 | 510 | 500 | 200 | 2 | 4 | 7 | VIII |

320.0 | 1.610 | 510 | 510 | 500 | 100 | 1 | 4 | 7 | VIII |

340.0 | 1.710 | 510 | 510 | 1000 | 100 | 1 | 5 | 7 | IX |

374.8 | 1.010 | 510 | 510 | 500 | 152 | 2 | 5 | 7 | IX |

[1/(m yr)] | [m/yr] | [m/yr] | [m/yr] | [m] | [m] | |||

= 1.3 = 0.5 | ||||||||

63.2 | 1.0010 | 510 | 5 | 0 | 100 | 1 | 0 | 0 |

94.9 | 1.5010 | 510 | 5 | 0 | 100 | 1 | 1 | 1 |

189.7 | 3.0010 | 510 | 5 | 0 | 100 | 1 | 3 | 1 |

389.6 | 6.1610 | 510 | 5 | 0 | 100 | 1 | 5 | 1 |

506.0 | 8.0010 | 510 | 5 | 0 | 100 | 1 | 5 | 1 |

594.5 | 9.4010 | 510 | 5 | 0 | 100 | 1 | 7 | 2 |

632.5 | 1.0010 | 510 | 5 | 0 | 100 | 1 | 7 | 2 |

948.7 | 1.5010 | 510 | 5 | 0 | 100 | 1 | 9 | 3 |

1264.9 | 2.0010 | 510 | 5 | 0 | 100 | 1 | 11 | 4 |

1897.4 | 3.0010 | 510 | 5 | 0 | 100 | 1 | 11 | 6 |

= 0.7 = 0.5 | ||||||||

10.3 | 2.0010 | 510 | 10 | 0 | 100 | 1 | 0 | 0 |

20.5 | 4.0010 | 510 | 10 | 0 | 100 | 1 | 1 | 1 |

30.8 | 6.0010 | 510 | 10 | 0 | 100 | 1 | 1 | 3 |

41.1 | 8.0010 | 510 | 10 | 0 | 100 | 1 | 1 | 3 |

43.7 | 8.5010 | 510 | 10 | 0 | 100 | 1 | 1 | 3 |

46.2 | 9.0010 | 510 | 10 | 0 | 100 | 1 | 2 | 3 |

51.4 | 1.0010 | 510 | 10 | 0 | 100 | 1 | 2 | 5 |

56.5 | 1.1010 | 510 | 10 | 0 | 100 | 1 | 2 | 5 |

= 1 = 0.7 | ||||||||

10.0 | 2.0010 | 510 | 5 | 0 | 100 | 1 | 0 | 0 |

50.2 | 1.0010 | 510 | 5 | 0 | 100 | 1 | 1 | 1 |

100.5 | 2.0010 | 510 | 5 | 0 | 100 | 1 | 1 | 3 |

125.6 | 2.5010 | 510 | 5 | 0 | 100 | 1 | 1 | 3 |

150.7 | 3.0010 | 510 | 5 | 0 | 100 | 1 | 2 | 3 |

175.8 | 3.5010 | 510 | 5 | 0 | 100 | 1 | 2 | 3 |

185.9 | 3.7010 | 510 | 5 | 0 | 100 | 1 | 2 | 5 |

= 1 = 0.3 | ||||||||

15.9 | 2.0010 | 510 | 10 | 0 | 100 | 1 | 0 | 0 |

47.8 | 6.0010 | 510 | 10 | 0 | 100 | 1 | 1 | 1 |

87.6 | 1.1010 | 510 | 10 | 0 | 100 | 1 | 1 | 3 |

119.4 | 1.5010 | 510 | 10 | 0 | 100 | 1 | 1 | 5 |

159.2 | 2.0010 | 510 | 10 | 0 | 100 | 1 | 1 | 7 |

199.1 | 2.5010 | 510 | 10 | 0 | 100 | 1 | 1 | 9 |

222.9 | 2.8010 | 510 | 10 | 0 | 100 | 1 | 1 | 9 |

254.8 | 3.2010 | 510 | 10 | 0 | 100 | 1 | 2 | 11 |

278.7 | 3.5010 | 510 | 10 | 0 | 100 | 1 | 2 | 11 |

318.5 | 4.0010 | 510 | 10 | 0 | 100 | 1 | 3 | 11 |