# Path selection in the growth of rivers

^{a}Lorenz Center, Department of Earth, Atmospheric and Planetary Sciences, Massachusetts Institute of Technology, Cambridge, MA 02139;^{b}Institut de Physique du Globe, 75252 Paris Cedex 05, France;^{c}Institute of Theoretical Physics, Faculty of Physics, University of Warsaw, 02-093, Warsaw, Poland

See allHide authors and affiliations

Edited by Grigory Isaakovich Barenblatt, University of California, Berkeley, CA, and approved September 9, 2015 (received for review July 21, 2014)

## Significance

The complex patterns of river networks evolve from interactions between growing streams. Here we show that the principle of local symmetry, a concept originating in fracture mechanics, explains the path followed by growing streams fed by groundwater. Although path selection does not by itself imply a rate of growth, we additionally show how local symmetry may be used to infer how rates of growth scale with water flux. Our methods are applicable to other problems of unstable pattern formation, such as the growth of hierarchical crack patterns and geologic fault networks, where dynamics remain poorly understood.

## Abstract

River networks exhibit a complex ramified structure that has inspired decades of studies. However, an understanding of the propagation of a single stream remains elusive. Here we invoke a criterion for path selection from fracture mechanics and apply it to the growth of streams in a diffusion field. We show that, as it cuts through the landscape, a stream maintains a symmetric groundwater flow around its tip. The local flow conditions therefore determine the growth of the drainage network. We use this principle to reconstruct the history of a network and to find a growth law associated with it. Our results show that the deterministic growth of a single channel based on its local environment can be used to characterize the structure of river networks.

As water flows, it erodes the land and produces a network of streams and tributaries (1⇓–3). Each stream continues to grow with the removal of more material, and evolves in a direction that corresponds to the water flux entering its head. The prediction of the trajectory of a growing channel and the speed of its growth are important for understanding the evolution of complex patterns of channel networks. Several models address their evolution and ramified structure. One, the Optimal Channel Networks model (4), is based on the concept of energy minimization and suggests a fractal network. The landscape evolution method and many diffusion-based models (5⇓–7) have also proven useful for modeling erosion and sediment transport. These models distinguish between two regimes: one is a diffusion-dominated regime where topographic perturbations are diminished, which leads to a smoother landscape and uniform symmetric drainage basins. In this case, the shape of a channel cannot deviate from a straight line. In the second regime, advection dominates, and channel incisions are amplified. The channel effectively continues to the next point that attracts the largest drainage basin, which corresponds to the direction where it receives the maximum water flux. These models nicely predict the formation of ridges and valleys and provide insight into the interaction between advective and diffusive processes (8, 9). However, they do not address directly the evolution of a single channel and do not explicitly address the nature of a growing stream based on its local environment.

Here we address two basic questions in the evolution and the dynamics of a growing channel: where it grows and at what velocity. We propose that, when streams are fed by groundwater, the direction of the growth of a stream is defined by the groundwater flow in the vicinity of the channel head. This theory is widely used in the framework of continuum fracture mechanics and accurately predicts crack patterns in different fracture modes, for both harmonic and biharmonic fields and for different stress singularities (10⇓–12). The theory, known as the principle of local symmetry, states that a crack propagates along the direction where the stress distribution is symmetric with respect to the crack direction (10, 13). We find an analog of this principle in the motion and growth of channels in a diffusive field. We argue that the trajectory of a stream is dictated by the symmetry of the field in the vicinity of its head. We also demonstrate how to determine a growth law that ties the water flux into a channel to an erosion process and the propagation velocity of a channel head.

## Principle of Local Symmetry

We study the case of channel growth driven by groundwater seepage as a process representative of channel formation and growth in a diffusing field (2, 14⇓⇓⇓⇓–19). The emergence of groundwater through the surface leads to erosion and the development of a drainage network (17). The flow of groundwater is described by Darcy’s law (2),*κ* the hydraulic conductivity, *p* the pressure in the fluid, *ρ* the fluid density, *z* the geometric height, and *g* the gravitational acceleration. By assuming only horizontal flow, the Dupuit approximation (20, 21) relates the water table height *P* and *κ*,*h* is a solution of the Poisson equation (21, 22). We assume that the hydraulic conductivity *κ* is constant. By rescaling the field, Eq. **2** becomes

Here we investigate how the groundwater flow controls the growth of a channel. In the vicinity of the channel head we can neglect the Poisson term in Eq. **3**, and the field can be approximated as (23, 24)*x* axis with boundary conditions*r* is the distance from the channel head and the channel is located at

Because the leading term in the expansion of Eq. **6** is symmetric with respect to *θ*, we do not expect it to influence the direction in which the channel grows. Thus, we must consider the subdominant term that breaks the symmetry and can therefore cause the stream to bend as it grows. Other terms in the expansion are negligible in the vicinity of the channel head.

In fracture mechanics, a crack maintains a symmetric elastic field around its tip to release the maximum stress as it propagates (10, 13). Inspired by this example, we suggest that a channel grows in the direction which maintains a locally symmetric groundwater flow. Accordingly, we formulate an analog of the principle of local symmetry as follows: A channel grows in the direction for which the coefficient

## Evaluation of the Principle of Local Symmetry

We now evaluate the principle of local symmetry (PLS) from three points of view. First, we justify its application to channel networks based on simple physical reasoning. We then show that the PLS is mathematically equivalent to assuming that a channel grows along the groundwater flow lines, an assumption that has been shown to be consistent with field observations (24). Finally, we describe a numerical procedure to grow a network according to this principle.

### Negative Feedback Induced by Groundwater.

As a channel grows, it moves the boundary of the network, and this changes the groundwater flow near the tip of the channel. Mathematically, the coefficients **6**] vary during the channel’s growth.

To understand how groundwater controls growth, we study the influence of the coefficient of the first two dominant terms (as *r* from the channel’s head, the amount of water collected through the right-hand side of the stream reads**7** and **8** indicate that more groundwater seeps into the stream through its right-hand bank (

### Equivalence with Geodesic Growth.

The geometry of a network draining a diffusion field depends on the dynamics of its growth. In particular, when two nascent tips grow off of a parent channel, the angle of this bifurcation depends on the growth rule. For instance, if a channel grows along the stream line that intersects its tip (Fig. 3*A*), it should bifurcate at

To grow numerically a network in a diffusion field, one typically needs to (*i*) solve the diffusion field around the network, (*ii*) evolve the network according to the local properties of the field near its tips, and repeat these operations (23). Whether the diffusion field satisfies the Poisson equation **3** or the Laplace equation **4** does not impact significantly the numerical procedure. However, complex analysis facilitates considerably the derivation of formal results about network growth in Laplacian fields.

In particular, the Laplacian field representing the groundwater flow defines a complex map from the physical plane to the upper half-plane (Fig. 3). This map encapsulates both the geometry of the network at a specific time, and the groundwater flow around it. Translating geodesic growth and local symmetry into this formalism, we find that they define the same growth rule (*SI Equivalence Between Geodesic Growth and Local Symmetry*).

This equivalence does not necessarily hold when the channel grows in a field which satisfies the Poisson equation. However, because the above demonstration is based on the local properties of the field, the source term of the Poisson equation comes as an external flux only. We therefore expect that the PLS is but a reformulation of geodesic growth, which accords with field observation. An important consequence of this reformulation is the ability to specify the growth direction in a precise, well-controlled manner, to which we now turn.

### Numerical Implementation.

We design a numerical method to calculate trajectories that explicitly maintain local symmetry, and compare its results to an analytic solution. We consider a simple case in which one channel grows in a confined rectangular geometry *SI Propagation of a Channel*.) We apply the following boundary conditions: a zero elevation at the bottom (

To validate our numerical implementation of the PLS, we compare our numerical trajectory to the evolution of a path in a Laplacian field according to the deterministic Loewner equation (26⇓–28). In the Loewner model, the properties of analytic functions in the complex plane are used to map the geometry into the complex half-plane or into radial geometry, and to find the solution for the field. Then, at each time step, a slit is added to the tip of the channel based on the gradient of the field entering the tip. Fig. 4 compares results from the two approaches. We find that the two solutions exhibit the same trajectory.

## SI Equivalence Between Geodesic Growth and Local Symmetry

Here we propose a detailed version of the reasoning briefly exposed in the main text. We first formulate the growth of a channel in a Laplacian field as the evolution of a conformal map, through the Loewner equation. We then establish the equivalence between geodesic growth and local symmetry in this framework.

### Conformal Map.

The diffusive field which surrounds the network satisfies Laplace’s equation. It is therefore the imaginary part of an analytical function **5**. By construction, the imaginary part of **4** at time *t*,

### Loewner Equation.

Following Loewner, we decompose the growth of a channel into a series of infinitesimal slit maps, each of which further increments the length of the channel (Fig. S1) (28).

Setting the length of each increment to zero, the continuous equation describing the growth of *N* channels in the physical plane reads (28)*t* for the *n*th channel in the mathematical plane. Note that the poles are the position of the tips in the mathematical plane. The growth factor sets the velocity at which a tip grows, whereas the motion of the pole sets the direction of its growth.

We define the inverse map **S3** with respect to *t* defines the evolution of the inverse map as*ω* and *t*, respectively.

The motion of the poles

### Geodesic Growth.

By definition, geodesic growth means that the channel’s tip at time *t*, i.e., *t*, and *φ*, Eq. **S5**, we find that geodesic growth implies

Expanding the Loewner function *φ* to second order near the pole [note that **S8**, we find, at leading order in

### PLS.

Local symmetry in the physical plane means that *g* up to the first order near the tip’s position *γ* becomes*g*, Eq. **S1**, the coefficient **6**, such that *f* up to the third order in the vicinity of the pole *p* becomes

Substitute the two expansions **S10** and **S11** in Eq. **S3**, and we find the relationship between the physical coefficients *c*. In particular, we find

Finally, applying the expansion **S11** into the Loewner equation **S4** for a single pole (*φ* controls the growth velocity. Eq. **S14** establishes that whenever a tip grows (**S9**).

## SI Propagation of a Channel

The evolution of a channel is defined by the field in its vicinity. To find the direction in which a channel tip evolves, we developed a numerical solver using Galerkin finite element discretization on a triangular grid (39), and solved the Poisson equation, Eq. **3**, or the Laplace equation, Eq. **4**, with the described boundary conditions. Then we iteratively add a small segment to a stream in different directions (Fig. S2). We obtain the field by using the numerical solver, and from the field of the vicinity of the studied stream, we find the coefficients of the expansion **6**. We accept the growth in the direction where the asymmetric coefficient becomes zero,

## Growth of a Real Stream Network

The evolution of a channel is defined by the field in the vicinity of the tip. However, this field is nonlocal and highly dependent on the boundary conditions imposed by the ramified network of the streams. In this section, the numerical method developed in the previous section is used to compute trajectories in more general settings where no analytic solutions are possible.

### Growth According to Local Symmetry.

We seek to determine if growth of a real stream network is consistent with the PLS. We study a network of seepage valleys located near Bristol, FL on the Florida Panhandle (17). The network is presented in Fig. 5. The valley network is obtained from a high-resolution LIDAR (Light Detection and Ranging) map with horizontal resolution of 1.2 m and vertical resolution of less than 5 cm (17). In this network, groundwater flows through unconsolidated sand above the relatively impermeable substratum, and into the streams (17, 29). Previous analyses indicate that the homogeneity of the sand is consistent with the assumption of constant *κ* (18, 19, 24, 30). The flow is determined by the Poisson equation **2**; thus, the network grows in a Poisson field (30).

We study the evolution of this network and check if the growth of the streams fulfills local symmetry. First, we set the boundary conditions; because the change in elevation along the Florida network is small (the median slope **3**, and find for each channel head the coefficient **6** that corresponds to the water flux entering the tip. Then, we remove a segment, *i*th tributary and propagate it forward to its original length in five small steps (to reduce numerical error). The growth of each stream is characterized by two variables: its growth rate and the direction of its growth. We assume that the velocity of a stream is proportional to the magnitude of the gradient of the field, raised to a power *η*:*i*th channel (*n* tributaries to be *n* meters. Each channel then grows in a direction that fulfills local symmetry, i.e., in the direction for which *β*, between the real trajectory of the stream and the reconstructed trajectory. We perform this calculation for 255 channel heads in the Florida network. We obtain a mean of *β* around zero with an SD of *η* give a similar distribution. We find that a mean around 0 of the angle error is consistent with a growth that fulfills local symmetry. However, some of the streams deviate significantly from their real growth direction, which may suggest that other factors account for their growth.

To evaluate the significance of the results, we suggest a null hypothesis in which the streams grow in the direction of the tangent regardless of the value of *φ*, between the tangent direction and the real trajectory, as shown in Fig. 6. We find that a growth according to the PLS reduces on average the error angle by *F* test (34) shows that the reduction of the variance using PLS is statistically significant, with *P*

### Growth Law.

To understand the deviations between the real and the calculated path, we hypothesize that the deviant streams grew in a different environment than currently exists, e.g., the neighboring tributaries were relatively undeveloped (or overdeveloped) when the studied stream reached its current location. To illustrate this idea, in Fig. 7 we show the trajectories of two streams with different velocities, and compare their evolution in Poisson field for different growth exponents *η*. One notices that for smaller *η* the slower streams are more likely to deviate from their real trajectory, but for higher *η* the faster streams change their course. Only when *η*), any errors will remain uncorrelated to the velocity of the streams.

Motivated by this reasoning, we study the correlation between the flux entering the tip, which we identify with *β* for different values of *η*. Retracting the network with different *η* creates different boundary conditions and influences the trajectory of the streams as they grow forward. For small *A*). Thus, as we grow the network forward, the deviation from the real trajectory will be larger for the slower streams, with small *β*, which suggests that

The importance of the growth exponent *η* is in the evolution of the network: Negative *η* will generate a stable network in which each perturbation, or small channel, will survive regardless of the water flux entering the tip. A positive *η* results in an unstable structure in which a small difference in the velocity of two competing channels is amplified and may lead to a screening mechanism and the survival only of the faster channel (28). Fig. 9 contains a schematic representation of this concept. The small positive exponent found for the stream network in Florida indicates that this network is unstable. This conclusion is consistent with the prediction of a highly ramified network.

## Summary

In summary, we offer a criterion for path selection of a stream in a diffusing field. We show that this criterion, which is based on the PLS (10, 13), predicts accurately the evolution of channels fed by groundwater. We suggest a method to infer the history of a real network by reconstructing it according to the PLS and evaluating errors for different growth laws. We parameterize the large-scale relationship between water flux and sediment transport with a single exponent and show that for the Florida network this growth exponent is roughly 0.7. We envision that our methods may also be applied to other problems, such as the growth of hierarchical crack patterns (35⇓–37) and geological fault networks (38), to provide a better understanding of their evolution.

## Acknowledgments

Inspiration for this work derived in part from visits to the Apalachicola Bluffs and Ravines Preserve. We thank D. Printiss and The Nature Conservancy for access to the Preserve. This work was supported by Department of Energy Grant FG02-99ER15004. P.S. acknowledges the support by the National Science Centre (Poland) under Research Grant 2012/07/E/ST3/01734.

## Footnotes

- ↵
^{1}To whom correspondence should be addressed. Email: ycohen{at}mit.edu.

Author contributions: Y.C., O.D., H.F.S., and D.H.R. designed research; Y.C., O.D., R.S.Y., P.S., and D.H.R. performed research; Y.C. and D.H.R. analyzed data; and Y.C., O.D., P.S., and D.H.R. wrote the paper.

The authors declare no conflict of interest.

This article is a PNAS Direct Submission.

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

## References

- ↵.
- Horton R E

- ↵.
- Dunne T

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

- ↵
- ↵
- ↵
- ↵
- ↵
- ↵Perron JT, Dietrich WE, Kirchner, JW (2008) Controls on the spacing of first-order valleys..
*J Geophys Res*113(F4):F04016 - ↵
- ↵
- ↵
- ↵
- ↵Dunne T (1969) Runoff production in a humid area. PhD thesis (Johns Hopkins University, Baltimore).
- ↵.
- Dunne T

- ↵Dietrich WE, Dunne, T (1993) The channel head..
*Channel Network Hydrology*, eds Beven K, Kirkby MJ (John Wiley and Sons, Chichester, UK), pp 175–219 - ↵Abrams DM, Lobkovsky AE, Petroff AP, Straub KM, McElroy B, Mohrig DC, Kudrolli A, Rothman DH (2009) Growth laws for channel networks incised by groundwater flow..
*Nat Geosci*2(3):193–196 - ↵
- ↵
- ↵.
- Bear J

- ↵.
- Dupuit J

- ↵.
- Polubarinova-Kochina PY

- ↵Petroff AP, Devauchelle O, Seybold H, Rothman DH (2013) Bifurcation dynamics of natural drainage networks..
*Philos Trans R Soc, A*371(2004):20120365 - ↵.
- Devauchelle O,
- Petroff AP,
- Seybold HF,
- Rothman DH

- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵.
- Howard AD,
- Kerby G

- ↵.
- Stock JD,
- Montgomery DR

*J Geophys Res: Solid Earth*104(B3):4983–4993 - ↵.
- Whipple KX,
- Tucker GE

*J Geophys Res: Solid Earth*104(B8):17661–17674 - ↵Press WH (2007).
*Numerical Recipes, 3rd Edition: The Art of Scientific Computing*(Cambridge Univ Press, Cambridge, UK) - ↵Skjeltorp AT, Meakin, P (1988) Fracture in microsphere monolayers studied by experiment and computer simulation..
*Nature*335(6189):424–426 - ↵
- ↵
- ↵Dolan JF, Bowman DD, Sammis CG (2007). Long-range and long-term fault interactions in southern california..
*Geology*35(9):855–858 - ↵.
- Kwon YW,
- Bang H

## Citation Manager Formats

## Article Classifications

- Physical Sciences
- Applied Physical Sciences