## New Research In

### Physical Sciences

### Social Sciences

#### Featured Portals

#### Articles by Topic

### Biological Sciences

#### Featured Portals

#### Articles by Topic

- Agricultural Sciences
- Anthropology
- Applied Biological Sciences
- Biochemistry
- Biophysics and Computational Biology
- Cell Biology
- Developmental Biology
- Ecology
- Environmental Sciences
- Evolution
- Genetics
- Immunology and Inflammation
- Medical Sciences
- Microbiology
- Neuroscience
- Pharmacology
- Physiology
- Plant Biology
- Population Biology
- Psychological and Cognitive Sciences
- Sustainability Science
- Systems Biology

# Channelization cascade in landscape evolution

Edited by Ignacio Rodriguez-Iturbe, Texas A&M University, College Station, TX, and approved December 6, 2019 (received for review July 11, 2019)

## Significance

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 describing the interplay between soil creep, runoff erosion, and tectonic uplift. We identify the critical conditions for the transition from a smooth to a channelized topography by means of a linear stability analysis and highlight striking similarities with fluid dynamic turbulence. The results shed light on the physical mechanisms responsible for observed landscape self-organization. The formation of regular prefractal networks reveals the tendency of the system to evolve toward optimal configurations typical of nonequilibrium complex systems.

## 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 partial differential equations coupling landscape evolution dynamics with a specific catchment area equation. By means of a linear stability analysis we identify the critical conditions triggering channel formation and the emergence of characteristic valley spacing. The ensuing channelization cascade, described by a dimensionless number accounting for diffusive soil creep, runoff erosion, and tectonic uplift, is reminiscent of the subsequent instabilities in fluid turbulence, while the structure of the simulated patterns is indicative of a 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 strongly influences the channelization cascade, underlying the nonlocal and nonlinear character of its dynamics.

The spatial distribution of ridges and valleys, including the formation of evenly spaced first-order valleys as well as more complex branching river networks (Fig. 1), is one of the most striking features of a landscape. It has long fascinated the scientific community, 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 (1⇓–3) and the coupled dynamics of water and sediment transport to identify stability conditions for incipient valley formation (4⇓–6), followed by the statistical description of river networks, including scaling laws and fractal properties of river basins (7⇓⇓–10), the related optimality principles (9, 11), and stochastic models (12⇓–14). 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 (13, 15, 16), 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 (17). To this purpose, landscape evolution models have been employed for the analysis of branching river networks (18, 19) 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 (20, 21) precluded the theoretical investigation of the underlying instabilities in relation to the leading geomorphological processes involved.

In this work, we focus on landscapes characterized by runoff erosion, expressed as a function of the specific drainage area a (21) 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 and allows us to describe the resulting ridge and valley patterns as a function of the relative proportions of diffusive soil creep, runoff erosion, and tectonic uplift. The nonlocal character of the equations makes the boundary conditions extremely important. On regular (i.e., square and rectangular) domains, simulations reveal a sequence of channel instabilities reminiscent of the laminar-to-turbulent transition (22⇓–24). The explicit mathematical structure makes it possible to perform a linear stability analysis of the coupled PDE system to identify the critical conditions for the first channel-forming instability. The subsequent branching sequence toward 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 networks of ridges and valleys, brought about by the regular boundary conditions, also reveals the tendency of the system to develop configurations suggestive of optimization principles (11) typical of nonequilibrium thermodynamics and complex systems (16, 25⇓⇓⇓⇓⇓⇓–32). Our analysis is different from recent interesting contributions on groundwater-dominated landscapes (33, 34), where branching and valley evolution is initiated at seepage points in the landscape.

## Landscape Evolution in Detachment-Limited Conditions

The time evolution of the surface elevation **1** becomes

The surface water flux q is linked to the continuity equation**3** can be simplified assuming steady-state conditions with constant, representative rainfall rate, **2** and **3** reduces to an equation for the specific catchment area a (21),

It is important to observe that the specific drainage area a has units of length and is related to the drainage area A as **5** instead of a, with several notable implications. The value of A is generally evaluated using numerical flow-routing algorithms [e.g., D8, D∞ (43)] which provide grid-dependent values of A. To correct for this, the drainage area A is often modified to account for the channel width (18, 41), but this results in approximations with arbitrary parameters. Conversely, the use of a avoids grid dependence of the resulting topography. Moreover, recasting the problem in terms of a consistent coupled system of PDEs makes it possible to analyze theoretically the landscape evolution process. As detailed below (*Materials and Methods*), an analytic solution for the steady-state hillslope profile can be derived (44) and then used as a basic state for a linear stability analysis to identify the critical conditions for the first channel formation and the characteristic valley spacing.

It is useful to nondimensionalize the system of Eqs. **4** and **5** to quantify the relative impact of soil creep, runoff erosion, and uplift on the landscape morphology. Using a typical length scale of the domain, l, and the parameters of Eqs. **4** and **5**, the following dimensionless quantities can be introduced: **5** becomes

## Results

### Organized Ridge and Valley Patterns.

Simulation results obtained by numerically solving Eqs. **4** and **5** over square domains with *Materials and Methods* for details) are shown in Fig. 2. The emerging ridge/valley patterns are classified in terms of Shreve order (used here as a measure of branching complexity) (3) and number of channels formed on each side of the domain. As can be seen from Eq. **7**, for *A*). As the runoff erosion coefficient is increased the system progressively develops 1, 3, and 5 channels on each side of the square domain for *B*–*D*). When *E*). As *F*–*I*), becoming of order 9 for the highest

As the resulting landscape changes from a smooth topography to a progressively more dissected one, the shape of the hypsometric curve varies from concave (i.e., slope decreases along the horizontal axis) to one with a convex portion for low elevations (Fig. 2*K*). 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 profile. As secondary branches develop, the hypsometric curve shifts to a concave/convex one, with the convex portion at lower altitudes becoming more evident as *K*).

The striking regularity of the drainage and ridge patterns induced by the simple geometry of the domain is reminiscent of regular prefractal structures [e.g., Peano basin (8, 9, 46⇓–48)] and is indicative of the fundamental role of boundary conditions due to the highly nonlocal control introduced by the drainage area term. The introduction of noise and irregular boundaries quickly breaks the regularity of the patterns (see results from numerical simulations obtained over progressively more irregular boundaries in *SI Appendix*, Fig. S10). The ridge and valley networks of Fig. 2 highly resemble figure 5 in ref. 31, where optimized tree-shaped flow paths were constructed to connect one point to many points uniformly distributed over an area. We further highlight similarities with the patterns obtained in ref. 30 by means of an erosion model where the global flow resistance is minimized.

### Effect of Runoff Erosion Laws.

The effect of different runoff erosion laws has been discussed in the literature (42) also in relation to climate, vegetation cover, and soil properties (49, 50). Their role was analyzed here by changing the values of the exponents n and m, as shown in Fig. 3.

When the value of n is different from unity, the resulting ridge/valley patterns depend on the uplift rate U, as in Eq. **7**. When n is increased, the system displays channelization and secondary branching for higher values of *A* and *B*), with a more dissected planar geometry characterized by narrower valleys and smaller junction angles (Fig. 3 *C*–*F*). A decrease in n leads to smoother geometries with wider valleys and the first secondary branching developing when only 3 channels per each side of the domain are present (Fig. 3 *G*–*J*). This results in a hypsometric curve with a more pronounced basal (i.e., at low altitudes) convexity for *SI Appendix*, Fig. S2).

As m is increased (Fig. 3 *K*–*N*) the system develops secondary branching when only 3 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 (*SI Appendix*, Fig. S2). Conversely, a decrease in m results in a more dissected landscape, with narrower valleys (Fig. 3 *O*–*R*) and a more pronounced transition of the hypsometric curve to a convex shape for low altitudes (*SI Appendix*, Fig. S2).

### Wide Rectangular Domains.

To assess boundary-condition effects on branching patterns we also considered very wide rectangular domains (*Materials and Methods*), around which we also performed a linear stability analysis. In our analogy with turbulent flows, the case of wide rectangular domains corresponds to the flow of viscous fluids between parallel plates (23, 24).

Results from the linear stability analysis are shown in Fig. 4. A critical value *C*), with the formation of progressively narrower valleys. Results from the linear stability analysis are in line with predictions from numerical experiments conducted over large rectangular domains, where the first channel instability occurs at *SI Appendix*, Fig. S9, as the water flow exponent m decreases, the value of

The numerical simulations confirm the results of the linear stability analysis and are in agreement with those of ref. 18. Fig. 5 compares the drainage patterns obtained as a function of *A*). After the first channelization, valleys tend to narrow as *B* and *C*). Further increasing the runoff erosion component provides progressively more dissected landscapes with the emergence of secondary branching (Fig. 5 *D*–*F*). As in turbulent flows larger Reynolds numbers produce smaller and smaller vortices; here increasing

The mean elevation profiles, computed as average elevation values along the x axis and neglecting the terminal parts of the domain to avoid boundary effects, are shown in Fig. 5 *G*–*L*. As the topography becomes progressively more dissected with increasing *G*–*L*). Such a behavior of the mean elevation profiles for increasing

The transition from a smooth to a channelized topography with increasing *Materials and Methods* for details). Fig. 5*P* shows the relationship between **11**) and deviate from it at *Materials and Methods* and figure 7.3 in ref. 53).

The effect of boundary conditions on the spatial regularity of ridge and valley patterns becomes especially apparent when comparing simulations with different aspect ratios. As can be seen in Fig. 5 *M*–*O*, when the domain size is slightly changed, the spatial organization of ridges and valleys is modified (e.g., the more regular pattern obtained for *SI Appendix*, Fig. S8). This suggests that some optimal domain length is needed to accommodate 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 Appendix*, Figs. S3–S7). This results in the formation of dislocation defects, as highlighted in the example of Fig. 5 *M*–*O*, as it is typical in nonlinear pattern-forming PDEs (52).

## 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 detachment-limited conditions. The sequence of instabilities is reminiscent of the subsequent bifurcations in fluid dynamic instabilities (23, 24, 52) and is captured by a dimensionless number (

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, with the characteristics of regular prefractals [e.g., the Peano basin (8, 46⇓–48)], is quickly lost as effects of noise and irregular boundaries are introduced (*SI Appendix*, Fig. S10). The shape of the hypsometric curve depends on the level of channelization and branching (54) 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. When diffusion dominates, hypsometric curves display a less pronounced basal convexity (54). A systematic analysis of real topographies in terms of statistics of hypsometry, branching angles, and characteristic spacing would help infer values of

It will also be interesting to explore the differences in transient dynamics between the hypsometries 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 presents a more pronounced basal convexity (2) regardless of the value of

## Materials and Methods

### Analytical Solutions for m = n = 1 .

To derive one-dimensional steady-state solutions of the coupled PDE system (Eqs. **4** and **5**) we consider a symmetric hillslope of length l in the x direction, with divide at *A*, *Inset*). Assuming a fixed elevation **4** and **5** for

### Linear Stability Analysis.

We studied the stability of the basic state (Eqs. **8** and **9**) to perturbations ã and *SI Appendix*, Fig. S9.

### Numerical Simulations.

Numerical simulations were performed using forward differences in time and centered difference approximations for the spatial derivatives, considering regular square grids of lateral dimension l, as well as rectangular domains with shape factor β, defined as the ratio between the domain dimensions in the y and x directions (i.e., *SI Appendix*, Fig. S1, were performed for different grid sizes to validate the independence of the resulting patterns on the grid resolution). Convex profiles were used as an initial condition. Over wide rectangular domains for

### Dimensional Analysis of the Channelization Transition.

We proceed similarly to the analysis of turbulence transition in pipes and channels. There the relationship between the friction factor ξ and the Reynolds number **10** as*P*, *Inset*).

### Data and Code Availability.

The 1-m resolution LiDAR data for Calhoun and Gabilan Mesa can be downloaded from the OpenTopography facility (https://opentopography.org). The code used for the numerical simulations is described in ref. 63 and available on GitHub (https://github.com/ShashankAnand1996/LEM).

## Acknowledgments

We acknowledge support from NSF Grants EAR-1331846 and EAR-1338694, and BP through the Carbon Mitigation Initiative at Princeton University. The useful comments of the anonymous reviewers are also gratefully acknowledged. LiDAR data for Calhoun and Gabilan Mesa were obtained from the National Center for Airborne Laser Mapping with support from the NSF (Grants EAR-1339015, EAR-1331846, and EAR-1043051) and retrieved from https://opentopography.org.

## Footnotes

- ↵
^{1}To whom correspondence may be addressed. Email: aporpora{at}princeton.edu.

Author contributions: S.B. and A.P. designed research; S.B., M.H., and C.C. performed research; S.B. and M.H. analyzed data; and S.B. and A.P. wrote the paper.

The authors declare no competing interest.

This article is a PNAS Direct Submission.

Data deposition: The 1-m resolution LiDAR data for Calhoun and Gabilan Mesa can be downloaded from the OpenTopography facility (https://opentopography.org). The code used for the numerical simulations is described in ref. 63 and available on GitHub (https://github.com/ShashankAnand1996/LEM).

This article contains supporting information online at https://www.pnas.org/lookup/suppl/doi:10.1073/pnas.1911817117/-/DCSupplemental.

Published under the PNAS license.

## References

- ↵
- R. E. Horton

- ↵
- A. N. Strahler

- ↵
- ↵
- T. R. Smith,
- F. P. Bretherton

- ↵
- D. S. Loewenherz

- ↵
- ↵
- D. G. Tarboton,
- R. L. Bras,
- I. Rodriguez-Iturbe

- ↵
- A. Marani,
- R. Rigon,
- A. Rinaldo

- ↵
- I. Rodríguez-Iturbe,
- A. Rinaldo

- ↵
- ↵
- R. Rigon,
- A. Rinaldo,
- I. Rodriguez-Iturbe,
- R. L. Bras,
- E. Ijjasz-Vasquez

- ↵
- ↵
- E. Somfai,
- L. M. Sander

- ↵
- R. Pastor-Satorras,
- D. H. Rothman

- ↵
- ↵
- ↵
- A. Fowler

- ↵
- J. T. Perron,
- W. E. Dietrich,
- J. W. Kirchner

- ↵
- ↵
- J. C. Gallant,
- M. F. Hutchinson

- ↵
- ↵
- S. B. Pope

- ↵
- P. G. Drazin,
- W. H. Reid

- ↵
- P. K. Kundu,
- I. M. Cohen,
- D. W. Dowling

- ↵
- ↵
- H. Ozawa,
- A. Ohmura,
- R. D. Lorenz,
- T. Pujol

- ↵
- L. M. Martyushev,
- V. D. Seleznev

- ↵
- ↵
- L. M. Sander,
- E. Somfai

- ↵
- M. R. Errera,
- A. Bejan

- ↵
- S. Lorente,
- W. Wechsatol,
- A. Bejan

- ↵
- A. Bejan

- ↵
- O. Devauchelle,
- A. P. Petroff,
- H. F. Seybold,
- D. H. Rothman

- ↵
- ↵
- W. E. Dietrich et al.

*Prediction in Geomorphology*(Geophysical Monograph Series, vol. 135, Blackwell Publishing Ltd, Oxford, UK, 2003), pp. 103–132. - ↵
- T. R. Smith

- ↵
- ↵
- ↵
- A. D. Howard

- ↵
- I. Rodríguez-Iturbe et al.

- ↵
- J. D. Pelletier

- ↵
- A. Chen,
- J. Darbon,
- J-M. Morel

- ↵
- D. G. Tarboton

- ↵
- S. Bonetti,
- D. D. Richter,
- A. Porporato

*Earth Surf. Process. Landforms*, https://doi.org/10.1002/esp.4694 (2019). - ↵
- A. Bejan,
- S. Lorente

- ↵
- B. B. Mandelbrot

- ↵
- ↵
- A. Flammini,
- F. Colaiori

- ↵
- D. R. Montgomery,
- G. Balco,
- S. D. Willett

- ↵
- L. E. L. Lowman,
- A. P. Barros

- ↵
- ↵
- ↵
- R. L. Panton

- ↵
- G. Willgoose,
- G. Hancock

- ↵
- S. Bonetti,
- A. Porporato

- ↵
- ↵
- M. Abramowitz,
- I. A. Stegun

- ↵
- C. G. Canuto,
- M. Y. Hussaini,
- A. Quarteroni,
- T. A. Zang

- ↵
- C. Camporeale,
- C. Canuto,
- L. Ridolfi

- ↵
- G. B. Chirico,
- A. W. Western,
- R. B. Grayson,
- G. Blöschl

- ↵
- K. E. Sweeney,
- J. J. Roering,
- C. Ellis

- ↵
- B. R. Munson,
- D. F. Young,
- T. H. Okiishi,
- W. W. Huebsch

- ↵
- S. K. Anand,
- M. Hooshyar,
- A. Porporato

## Citation Manager Formats

## Sign up for Article Alerts

## Article Classifications

- Physical Sciences
- Environmental Sciences

## Jump to section

## You May Also be Interested in

*Ikaria wariootia*represents one of the oldest organisms with anterior and posterior differentiation.