# River landscapes and optimal channel networks

Contributed by Andrea Rinaldo, May 10, 2018 (sent for review March 14, 2018; reviewed by Armaghan Abed-Elmdoust, Angelika Steger, and Ramarathnam Venkatesan)

## Significance

Optimal channel networks (OCNs) are a well-studied static model of river network structures. We present exact results showing that every OCN is a natural river tree where a landscape exist such that the flow directions are always directed along its steepest descent. We characterize the family of natural river trees in terms of certain forbidden structures, called “$k$-path obstacles.” We thus determine conditions for river landscapes to imply a rooted tree as its network of topographic gradients. Our results are significant in particular for applications where OCNs may be used to produce statistically identical replicas of realistic matrices for ecological interactions.

## Abstract

We study tree structures termed optimal channel networks (OCNs) that minimize the total gravitational energy loss in the system, an exact property of steady-state landscape configurations that prove dynamically accessible and strikingly similar to natural forms. Here, we show that every OCN is a so-called natural river tree, in the sense that there exists a height function such that the flow directions are always directed along steepest descent. We also study the natural river trees in an arbitrary graph in terms of forbidden substructures, which we call $k$-path obstacles, and OCNs on a $d$-dimensional lattice, improving earlier results by determining the minimum energy up to a constant factor for every $d\ge 2$. Results extend our capabilities in environmental statistical mechanics.

### Sign up for PNAS alerts.

Get alerts for new articles, or get an alert when an article is cited.

River networks can be viewed as rooted trees that, when extracted from fluvial landscapes by topographic steepest descent directions, show deep similarities of their parts and the whole, often across several orders of magnitude, despite great diversities in their geology, exposed lithology, vegetation, and climate (1). The large related body of observational data provides quintessential examples of real physical phenomena that can be effectively modeled by using graph theory (2, 3). River networks are in fact the loopless patterns formed by fluvial erosion over a drainage basin. Such patterns, referred to as “spanning trees,” are such that a directed graph route exists for flows from every site of the catchment to an outlet and are strictly related to the topographical surface whose gradients determine the edges (the drainage directions). A specific class of spanning trees, called optimal channel networks (OCNs), was obtained by minimizing a specific functional (4, 5), later shown to be an exact property of the stationary solutions of the general equation describing landscape evolution (6, 7). The static properties and the dynamic origins of the scale-invariant structures of OCNs proved remarkable (1). OCNs are suboptimal (that is, dynamically accessible given initial conditions and quenched randomness frustrating the optimum search) configurations of a spanning network mimicking landscape evolution and network selection (

*SI Appendix*). Empirical and theoretical works have generally focused on the two-dimensional (2D) case, although recently (inspired by vascular systems), three-dimensional (3D) settings have also been analyzed for an OCN embedded in a lattice (8). Several exact results have been derived for OCNs (6–15).A model for a river network is obtained by taking a reasonably dense set of points on a terrain and then joining each point to a nearby point downhill. Experimental observations suggest that OCNs should satisfy the following properties: (where $N\left(v\right)=\left\{u:uv\in E\left(G\right)\right\}$ denotes the set of neighbors of node $v$. It is easy to see that the set of edges along which water flows is acyclic; in a single river basin, it will be an oriented spanning tree, with all water flowing to a unique root. Note that the notation $N\left(v\right)\cup \left\{v\right\}$ also includes the possibility of finding nodes with no outlet. The second empirical property provides a connection between the gradient of the descent and the landscape-forming flux of water flowing along the path. More precisely, in the (geological) steady state, if the landscape-forming water flow rate ${m}_{v}$ at a site $v$ increases, then the heights of the sites are modified by water erosion in such a way as to keep the product ${m}_{v}^{1/2}\mathrm{\Delta}{h}_{v}$ constant. The exponent $1/2$ is experimentally determined by local slope-area plots where one assumes that landscape-forming discharges are proportional to total contributing area and refers to the fluvial domain, thereby neglecting unchanneled parts of the landscape (1). Such an exponent also emerges as the leading term of a general nonlocal erosion term under reparametrization invariance and a small-gradient approximation (7). Thus, the first rule microscopically ensures that no cycles are created, while the second rule introduces a feedback between river structure and landscape erosion (1, 16–18).

*i*) each portion of the drainage basin has a single output; and (*ii*) a geomorphological stability condition holds. OCNs are defined on a finite graph $G$, consisting of a set $V\left(G\right)$ of vertices (or nodes), which correspond to sites in a drainage river basin, and a set $E\left(G\right)$ of edges between nearby sites. Moreover, each node $v$ is endowed with a height, ${h}_{v}$. Water from each site/node $v$ is entirely directed along the steepest descent—that is, toward the neighbor $u$ of $v$ which maximizes ${h}_{v}-{h}_{u}$—and we write$$\mathrm{\Delta}{h}_{v}=\underset{u\in N\left(v\right)\cup \left\{v\right\}}{max}\left\{{h}_{v}-{h}_{u}\right\},$$

Here, we present analytical results about the properties of the structure that minimizes the total gravitational energy loss. In particular, we will give a theoretical characterization of river networks in terms of forbidden substructures, which provides a simple way of determining which of the spanning trees of a graph are natural river trees. We will also study OCNs embedded in a $d$-dimensional lattice, proving, for each dimension $d\ge 2$, upper and lower bounds on the energy that differ by only a constant factor. In fact, river trees are 2D, while their landscapes are 3D objects where a drop in elevation is assigned along edges as a function of connectivity. $d$-dimensional optimal constructions are of notable interest, like, e.g., spanning vascular networks ($d=3$) delivering metabolites to a given volume (13). In a vascular system of size $L$, metabolic rates scale like ${L}^{d}$, while mass $M$ scales as $M\propto {L}^{d+1}$, leading to metabolism scaling as ${M}^{d/\left(d+1\right)}$ (13). Constructions in $d>3$ may consider time as an extra dimension [see, e.g., the work on Abelian sandpile (19) models of the Bethe lattice (20)]. Great interest thus exists for OCNs on general graphs for which important applications may arise in the future. We thus begin with a formal description of OCNs, and then we present the analytical results concerning the characterization of river trees and the energy of an OCN embedded in lattices of any dimension $d$.

For each constant $\gamma \in \left(\mathrm{0,1}\right)$, we can associate to the (rooted) spanning trees of a graph $G$ an energy function of the following form:where ${A}_{v}\left(T\right)$ is the number of vertices $u\in V\left(G\right)$ such that the path in $T$ from $u$ to the root contains the vertex $v$ (in other words, the number of vertices in the subtree of $T$ “rooted” at $v$). A spanning tree of $G$ which minimizes the energy in Eq.

$${\mathcal{E}}_{\gamma}\left(T\right)\u2254\sum _{v\in V\left(G\right)}{A}_{v}{\left(T\right)}^{\gamma},$$

[1]

**1**is called an OCN of $G$.In the case of river networks, if we imagine an open system where injection occurs at each site (vertex) at rate $1$, and flowing along the edges of $T$ until it reaches the root, then ${A}_{v}\left(T\right)$ represents the total flow of water out of $v$. The energy function defined in Eq. where ${m}_{v}$ is the mass of water leaving vertex $v$, so ${m}_{v}$ is proportional to ${A}_{v}\left(T\right)$. Furthermore, experimental evidence suggests the following relationship between ${A}_{v}\left(T\right)$ and the difference of heights of two adjacent nodes: $\mathrm{\Delta}{h}_{v}\propto {A}_{v}{\left(T\right)}^{-1/2}$ (1). We therefore obtain:Replacing the exponent $1/2$ with the parameter $\gamma \in \left(\mathrm{0,1}\right)$, and ignoring the (unimportant) constant factor, gives Eq.

**1**has its origins in the exact result that steady-state configurations of river networks minimize the total loss of gravitational energy. This corresponds to:$${\mathcal{E}}_{\gamma}\left(T\right)=\sum _{v\in V\left(G\right)}{m}_{v}g\mathrm{\Delta}{h}_{v},$$

[2]

$${\mathcal{E}}_{\gamma}\left(T\right)\propto \sum _{v\in V\left(G\right)}{A}_{v}{\left(T\right)}^{1/2}.$$

[3]

**1**. When $\gamma =1$ the class of directed networks that minimize the mean distance to the outlet is obtained, whose energy function is the same (6, 16). This is true if the functional in Eq.**1**is replaced by a more general one where $T$ is substituted with the graph $G$ itself, including all edges. For such functional, every tree is a local minimum (7), and minimization is carried out under the constraint that locally (for each node) and globally a conservation law holds matching injection rates and inflows/outflows (6, 7, 16).Given a rooted graph $G$ (with one or more roots), we say that $T$ is a “rooted spanning forest” of $G$ if $T$ is an acyclic subgraph of $G$, with one root of $G$ in each component of $T$, and all edges oriented toward the corresponding root. The following natural definition will play an important role in our investigation of the OCN(s) of a graph.

## Definition.

Let $G$ be a rooted graph, and $T$ be a rooted spanning forest of $G$. We say that $T$ is a “natural river spanning forest” of $G$ if there is a height function $f:V\left(G\right)\to \mathbb{R}$ such that if $v\in V\left(G\right)$ is a nonroot, and $u\in N\left(v\right)$, then $v\to u$ is an edge of $T$ if and only if

$$f\left(u\right)=min\{f\left(x\right):x\in N\left(v\right)\cup \left\{v\right\}\}.$$

[4]

We will assume throughout that the values of $f$ are distinct, and hence that the minimum in Eq.

**4**occurs at a unique vertex $x$. We remark that we will usually be interested in the case when $G$ has exactly one root, in which case we say that $T$ is a “natural river spanning tree” of $G$.Note that (

*i*) every nonroot vertex of $T$ has out-degree exactly $1$, (*ii*) $f\left(v\right)$ strictly decreases along every directed path in $T$, and (*iii*) the function $f$ attains its minimum at one of the roots. In*Theorem 1*, we will provide a simple characterization of the natural river spanning trees of a graph in terms of forbidden substructures. We will also show (*Theorem 3*) that every OCN in a graph is also a natural river spanning tree of that graph. As a consequence, to study (either analytically or numerically) the tree structures that minimize the energy function ${\mathcal{E}}_{\gamma}\left(T\right)$ defined in Eq.**1**, it suffices to consider natural river trees.## Which Networks Are River Trees?

The definition of a natural river spanning tree creates various constraints on its possible shape and structure; in this section we shall consider these constraints.

## Definition.

Let $G$ be a rooted graph, and let $T$ a spanning tree of $G$. A $k$-path obstacle for $T$ consists of a sequence of directed paths $\left({\eta}_{1},\dots ,{\eta}_{k}\right)$ in $T$, with the following property: For each $i{\in \mathbb{Z}}_{k}$, the last vertex of ${\eta}_{i}$ is connected to the first vertex of ${\eta}_{i+1}$ by an edge of $G$ that is not in $T$.

In other words, if ${\eta}_{i}$ is the path ${a}_{i}\to {c}_{i}\left(1\right)\to \cdots \to {c}_{i}\left({\ell}_{i}\right)\to {b}_{i}$ in $T$, andthen we say that $T$ contains a $k$-path obstacle (Fig. 1). We emphasize that the paths ${\eta}_{i}$ do not need to be vertex-disjoint (and it will be convenient in the proofs below to allow them to intersect); however, if $k$ is minimal such that $T$ contains a $k$-path obstacle, then the paths in any $k$-path obstacle will be disjoint. Observe that we will always consider the indices modulo $k$, so that (for example) ${b}_{0}={b}_{k}$ and ${b}_{1}={b}_{k+1}$, and that if $k=1$ then necessarily ${\ell}_{1}>0$.

$${b}_{1}{a}_{2},{b}_{2}{a}_{3},\dots ,{b}_{k}{a}_{1}\in E\left(G\right)\backslash E\left(T\right),$$

Fig. 1.

It is not difficult to show that a natural river tree cannot contain a $k$-path obstacle for any $k\in \mathbb{N}$; our first theorem states that this simple necessary condition is also sufficient, and hence characterizes the natural river trees.

## Theorem 1.

*Let*$G$

*be a finite graph with a single root. A spanning tree*$T$

*of*$G$

*is a natural river spanning tree of*$G$

*if, and only if, it has no path obstacle.*

Looking at Fig. 1, it may be deduced that a natural river tree by definition cannot contain any three-path obstacle (i.e., cannot contain the red dashed edges). The same reasoning can be extended to any natural river tree and value of $k$. We will give here only a formal proof of the necessary and sufficient conditions for $T$ to be a natural river network.

## Proof.

We will show first that, for each $k\ge 1$, a natural river tree $T$ cannot contain a $k$-path obstacle. Indeed, let $f:V\left(G\right)\to \mathbb{R}$ be the height function associated to $T$, and suppose that there exist paths ${\eta}_{1},\dots ,{\eta}_{k}$ in $T$, where ${\eta}_{i}$ is the path ${a}_{i}\to {c}_{i}\left(1\right)\to \cdots \to {c}_{i}\left({\ell}_{i}\right)\to {b}_{i}$ in $T$, and ${b}_{1}{a}_{2},{b}_{2}{a}_{3},\dots ,{b}_{k}{a}_{1}\in E\left(G\right)\backslash E\left(T\right)$. Since ${a}_{i}\to {c}_{i}\left(1\right)$ (or ${a}_{i}\to {b}_{i}$ if ${\ell}_{i}=0$) is in $T$, and ${a}_{i}{b}_{i-1}\in E\left(G\right)$, it follows from the definition that $f\left({b}_{i}\right)\le f\left({c}_{i}\left(1\right)\right)<f\left({b}_{i-1}\right)$ for each $i{\in \mathbb{Z}}_{k}$, where ${b}_{0}\u2254{b}_{k}$. But then,and this contradiction proves the claim.

$$f\left({b}_{1}\right)>f\left({b}_{2}\right)>\cdots >f\left({b}_{k}\right)>f\left({b}_{1}\right),$$

[5]

We now turn to our main task: showing that if $T$ has no path obstacle, then it is a natural river tree. Given such a rooted spanning tree $T$ in $G$, we introduce a directed graph $D$ with vertex set $V\left(G\right)$ and

$$\begin{array}{ccc}\hfill E\left(D\right)=\{x\to y:\text{either}x\to y\hspace{0.17em}\text{is an edge of}\hspace{0.17em}T\text{, or}& \hfill & \hfill \\ \hfill \exists z\hspace{0.17em}\text{such that}\hspace{0.17em}xz\in E\left(G\right)\backslash E\left(T\right)\text{and}\left(z\to y\right)\in E\left(T\right)\}.& \hfill & \hfill \end{array}$$

## Claim.

*The directed graph*$D$

*is acyclic.*

## Proof.

Suppose that $D$ contains a directed cycle with exactly $k$ edges that are not edges of $T$, and note that $k\ge 1$, since $T$ is a tree. Let us denote the vertices of the cycle as follows:where ${b}_{1},\dots ,{b}_{k}$ are the vertices whose out-neighbor in the cycle is not an edge of $T$, and ${\ell}_{1},\dots ,{\ell}_{k}\ge 0$. Thus, there exist vertices ${a}_{1},\dots ,{a}_{k}$ such that, for each $i{\in \mathbb{Z}}_{k}$, we have ${b}_{i-1}{a}_{i}\in E\left(G\right)\backslash E\left(T\right)$ (where ${b}_{0}\u2254{b}_{k}$), and ${a}_{i}\to {c}_{i}\left(1\right)$ (or ${a}_{i}\to {b}_{i}$ if ${\ell}_{i}=0$) is an edge of $T$. This gives a $k$-path obstacle, contradicting our assumption.$\square $

$$\begin{array}{ccc}\hfill {c}_{1}\left(1\right)\to \cdots \hspace{0.17em}\to {c}_{1}\left({\ell}_{1}\right)\to {b}_{1}\to {c}_{2}\left(1\right)\to \cdots \to {c}_{2}\left({\ell}_{2}\right)\to & \hfill & \hfill \\ \hfill {b}_{2}\to \cdots \to {b}_{k-1}\to {c}_{k}\left(1\right)\to \cdots \to {c}_{k}\left({\ell}_{k}\right)\to {b}_{k}\to {c}_{1}\left(1\right),& \hfill & \hfill \end{array}$$

It follows from the claim that the directed graph $D$ can be extended to a linear ordering $<$ on $V\left(G\right)$ (i.e., a linear ordering with $x>y$ for every edge $x\to y$ of $D$). Choose a height function $f$ such that $f\left(x\right)<f\left(y\right)$ if and only if $x<y$; we claim that $T$ is a natural river spanning tree with this height function.

Indeed, let $v\in V\left(G\right)$ be a nonroot, and note that there is a unique $u\in N\left(v\right)$ such that $v\to u$ is an edge of $T$, since all edges of $T$ are directed toward the root. Note that $v\to u$ is an edge of $D$, by definition, so $f\left(u\right)<f\left(v\right)$. Now let $u\hspace{0.17em}\ne \hspace{0.17em}w\in N\left(v\right)$, and observe that $w\to u$ is an edge of $D$ (since $wv\in E\left(G\right)\backslash E\left(T\right)$ and $v\to u$ is an edge of $T$), so $f\left(u\right)<f\left(w\right)$. Hence,as required. It follows that a spanning tree of $G$ with no path obstacle is a natural river spanning tree of $G$, completing the proof of

$$f\left(u\right)=min\{f\left(x\right):x\in N\left(v\right)\cup \left\{v\right\}\},$$

*Theorem 1*.$\square $Recall that an OCN (with parameter $\gamma $) in a graph $G$ is a spanning tree of $G$ that minimizes the energy ${\mathcal{E}}_{\gamma}\left(T\right)$. Our next theorem states that such a tree does not contain a path obstacle.

## Theorem 2.

*Let*$G$

*be a finite graph with a single root, let*$T$

*be a spanning tree of*$G$

*, and let*$k\in \mathbb{N}$

*and*$\gamma \in \left(\mathrm{0,1}\right)$

*. If*$T$

*is an OCN in*$G$

*with parameter*$\gamma $

*, then*$T$

*does not contain a*$k$

*-path obstacle.*

The proof of this theorem is presented in

*Materials and Methods*.*Theorems 1*and*2*immediately imply the following fundamental fact.## Theorem 3.

*Every OCN is a natural river spanning tree.*

## Energy of OCNs in $d$-Dimensional Grids.

Here, we will study in more detail the energy of an OCN in a $d$-dimensional lattice. In particular, we will focus on the case of the $d$-dimensional grid ${\left[n\right]}^{d}$, with root $\left(1,\dots ,1\right)$, and edges between vertices at ${\ell}_{\infty}$-distance $1$, although it will be clear from the proofs below that our results can be easily extended to other lattices. For each $n,d\in \mathbb{N}$ and $\gamma \in \left(\mathrm{0,1}\right)$, definethe energy of an OCN in the $d$-dimensional grid. The following theorem determines $T\left(n,d,\gamma \right)$ up to a constant factor (depending on $d$ and $\gamma $) for all values of the parameters. It extends results of Colaiori et al. (12), who proved the lower bounds for $d=2$ and $\gamma \hspace{0.17em}\ne \hspace{0.17em}1/2$ (in fact, on an orthogonal grid), and improves a result of Briggs and Krishnamoorthy (8), who obtained the lower bound $3{n}^{2}/2-O\left(n\right)$ in the case $d=2$ and $\gamma =1/2$.

$$T\left(n,d,\gamma \right)\u2254min\{{}_{\gamma}\left(T\right):T\mathrm{\text{is a spanning tree of}}{\left[n\right]}^{d}\},$$

## Theorem 4.

*For every integer*$d\ge 2$

*, we have*

$$T\left(n,d,\gamma \right)=\left\{\begin{array}{cc}\mathrm{\Theta}\left({n}^{d}\right),\hfill & if\hspace{0.17em}0<\mathrm{\gamma}<1-1/d,\hfill \\ \mathrm{\Theta}\left({n}^{d}\mathrm{log}n\right),\hfill & if\hspace{0.17em}\mathrm{\gamma}=1-1/d,\hfill \\ \mathrm{\Theta}\left({n}^{1+d\gamma}\right),\hfill & if\hspace{0.17em}1-1/d<\mathrm{\gamma}<1,\hfill \end{array}\right.$$

[6]

*as*$n\to \infty $.

## Proof.

Let us first prove the upper bounds, which follow from a straightforward recursive construction. We may assume for simplicity that $n={2}^{k}$, and consider ${2}^{d}$ copies of the construction ${T}_{n/2}$ for ${\left[n/2\right]}^{d}$, oriented so that ${2}^{d}-1$ of them have their root near the center of ${\left[n\right]}^{d}$, and the last one has its root at $\left(1,\dots ,1\right)$, Fig. 2. Join the vertex $\left(n/2,\dots ,n/2\right)$ in this last subcube to each of the roots of the other subcubes, which are each of the form $\left(n/2,\dots ,n/2\right)+\mathbf{x}$ for some $\mathbf{x}\in {\left\{\mathrm{0,1}\right\}}^{d}$, so as to form a spanning tree ${T}_{n}$ in ${\left[n\right]}^{d}$.

Fig. 2.

We claim that the subtree ${T}_{n}$ of ${\left[n\right]}^{d}$ constructed in this way contains the diagonal edge $\left(x+1,\dots ,x+1\right)\to \left(x,\dots ,x\right)$ for each $1\le x<n$. Indeed, this holds (trivially) when $n=1$, and by induction it follows for all $n={2}^{k}$, since we add the diagonal edge $\left(n/2+1,\dots ,n/2+1\right)\to \left(n/2,\dots ,n/2\right)$ in the induction step, and all other diagonal edges are in ${T}_{n}$ by the induction hypothesis. It follows that all of the mass from ${2}^{d}-1$ copies of ${\left[n/2\right]}^{d}$ flows down the diagonal of the final copy, and hence we have the following recursive formula:Since ${\left(y+z\right)}^{\gamma}\le {y}^{\gamma}+{z}^{\gamma}$, and noting that $\left(n/2\right){\left({n}^{d}-{\left(n/2\right)}^{d}\right.}^{\gamma}\le {n}^{1+d\gamma}$, it follows thatand, hence, noting that ${\mathcal{E}}_{\gamma}\left({T}_{1}\right)=1$, and recalling that $n={2}^{k}$, we obtainThis inequality immediately implies the claimed upper bounds. Indeed, if $1+d\gamma <d$ then the right-hand side of Eq. If $1+d\gamma >d$ then it is a decreasing geometric series, soFinally, if $1+d\gamma =d$ then all terms are equal, soas required.

$$\begin{array}{ccc}\hfill {\mathcal{E}}_{\gamma}\left({T}_{n}\right)& =\hfill & {2}^{d}\cdot {\mathcal{E}}_{\gamma}\left({T}_{n/2}\right)+\hfill \\ \hfill & \hfill & +\sum _{x=1}^{n/2}{\left(\right({A}_{\left(x,\dots ,x\right)}\left({T}_{n/2}\right)+{n}^{d}-{\left(n/2\right)}^{d}}^{\gamma}+\hfill \\ \hfill & \hfill & -{A}_{\left(x,\dots ,x\right)}{\left({T}_{n/2}\right)}^{\gamma}).\hfill \end{array}$$

[7]

$${\mathcal{E}}_{\gamma}\left({T}_{n}\right)\le {2}^{d}\cdot \hspace{0.17em}{\mathcal{E}}_{\gamma}\left({T}_{n/2}\right)+{n}^{1+d\gamma},$$

$$T\left(n,d,\gamma \right)\hspace{0.17em}\le \hspace{0.17em}{\mathcal{E}}_{\gamma}\left({T}_{n}\right)\hspace{0.17em}\le \hspace{0.17em}{n}^{1+d\gamma}+\cdots +{2}^{\ell d}{\left(\frac{n}{{2}^{\ell}}\right.}^{1+d\gamma}+\cdots +{n}^{d}.$$

[8]

**8**is an increasing geometric series, so$$T\left(n,d,\gamma \right)\u2a7d\frac{{n}^{d}}{1-{2}^{1+d\gamma -d}}=O\left({n}^{d}\right).$$

$$T\left(n,d,\gamma \right)\u2a7d\frac{{n}^{1+d\gamma}}{1-{2}^{d-1-d\gamma}}=O\left({n}^{1+d\gamma}\right).$$

$$T\left(n,d,\gamma \right)\u2a7d{n}^{d}{\mathrm{log}}_{2}\hspace{0.17em}n+1=O\left({n}^{d}\mathrm{log}n\right),$$

We next turn to the lower bounds. Observe first that the lower boundholds trivially for every spanning tree $T$ of ${\left[n\right]}^{d}$, and every $\gamma \in \left(\mathrm{0,1}\right)$, since we have ${A}_{v}\left(T\right)\ge 1$ for every vertex $v$. Hence, $T\left(n,d,\gamma \right)\ge {n}^{d}$ holds always.

$${\mathcal{E}}_{\gamma}\left(T\right)=\sum _{v\in V\left({\left[n\right]}^{d}\right)}{A}_{v}{\left(T\right)}^{\gamma}\u2a7e\sum _{v\in V\left({\left[n\right]}^{d}\right)}1={n}^{d},$$

For the other lower bounds, we will find it useful to work in a more general setting, assuming for simplicity that $n={3}^{k}$. Let ${\left[n\right]}_{\u25cb}^{d}$ denote a copy of the $d$-dimensional grid in which every vertex on the boundary is a root and definewhere each component of $F$ is assumed to have a single root (which lies on the boundary of ${\left[n\right]}^{d}$), and the energy ${\mathcal{E}}_{\gamma}\left(F\right)$ is defined to be the sum of the energies of the components. Note that $F\left(n,d,\gamma \right)\le T\left(n,d,\gamma \right)$, since a spanning tree in ${\left[n\right]}^{d}$ contains a rooted spanning forest of ${\left[n\right]}_{\u25cb}^{d}$. We claim that

$$\begin{array}{llll}\hfill & F\left(n,d,\gamma \right)\u2254min{\left\{\mathcal{E}\right.}_{\gamma}\left(F\right):F\mathrm{\text{is a rooted}}\hfill & \hfill & \hfill \\ \hfill & \mathrm{\text{spanning forest of}}{\left[n\right]}_{\u25cb}^{d}\},\hfill & \hfill & \hfill \end{array}$$

$$F\left(3n,d,\gamma \right)\ge {3}^{d}\cdot F\left(n,d,\gamma \right)+\hspace{0.17em}\frac{\gamma}{{3}^{d\left(1-\gamma \right)}}\cdot {n}^{1+d\gamma}.$$

[9]

To prove Eq. for any $x>0$. Therefore, moving the mass of the center subcube to the boundary increases the energy by at leastand, hence, we obtain Eq.

**9**, we partition ${\left[3n\right]}^{d}$ into ${3}^{d}$ subcubes of size ${\left[n\right]}^{d}$. In each subcube we must send the mass to the boundary of this subcube, and this gives a contribution of at least ${3}^{d}\cdot F\left(n,d,\gamma \right)$ to the energy. Additionally, the middle subcube has to send its mass from its boundary to the boundary of the large cube ${\left[3n\right]}^{d}$, and therefore each drop of water must pass through at least $n$ additional vertices on its way to a root. Note that, since ${A}_{v}\left(T\right)\le {\left(3n\right)}^{d}$ and $\gamma <1$, we have$${A}_{v}{\left(T\right)}^{\gamma}-{\left({A}_{v}\left(T\right)-x\right.}^{\gamma}\ge \gamma \cdot {A}_{v}{\left(T\right)}^{\gamma -1}\cdot x\ge \gamma \cdot {\left(3n\right)}^{d\gamma -d}\cdot x,$$

$${n}^{d}\cdot n\cdot \gamma \cdot {\left(3n\right)}^{d\gamma -d}=\frac{\gamma}{{3}^{d\left(1-\gamma \right)}}\cdot {n}^{1+d\gamma},$$

**9**(Fig. 3).Fig. 3.

Setting $C=\gamma \cdot {3}^{-\left(d+1\right)}$, and noting that $F\left(1,d,\gamma \right)=1$, we obtainwhere in the penultimate term $\ell =k-2$. If $1+d\gamma <d$ then this gives only a minor improvement over the trivial lower bound ${n}^{d}$, but for $1+d\gamma >d$ the improvement is more substantial. Indeed, it implies thatas $n\to \infty $. Finally, if $1+d\gamma =d$, then we obtainas required.$\square $

$$\begin{array}{ccc}\hfill T\left(n,d,\gamma \right)\ge F\left(n,d,\gamma \right)\ge C\cdot {n}^{1+d\gamma}+\cdots & \hfill & \hfill \\ \hfill \cdots +C\cdot {3}^{\ell d}{\left(\frac{n}{{3}^{\ell}}\right.}^{1+d\gamma}+\cdots +{n}^{d},& \hfill & \hfill \end{array}$$

$$T\left(n,d,\gamma \right)\u2a7e\frac{\gamma +o\left(1\right)}{{3}^{d+1}}\cdot \frac{{n}^{1+d\gamma}}{1-{3}^{d-1-d\gamma}}=\mathrm{\Omega}{n}^{1+d\gamma},$$

$$T\left(n,d,\gamma \right)\u2a7e\frac{\gamma}{{3}^{d+1}}\cdot {n}^{d}{\mathrm{log}}_{3}\hspace{0.17em}n-2=\mathrm{\Omega}{n}^{d}\mathrm{log}n,$$

An interesting (and likely difficult) challenge would be to decrease the implicit constant factors between the upper and lower bounds in

*Theorem 4*to factors of $1+o\left(1\right)$, and hence to determine $T\left(n,d,\gamma \right)$ asymptotically.## Discussion

When one considers only the structure of the network (without imposing a height function on the vertices), the requirement that a river network at stationarity minimizes total energy dissipation takes a simple and elegant form: We need to minimize the function

$${\mathcal{E}}_{\gamma}\left(T\right)\propto \sum _{v\in V\left(G\right)}{A}_{v}{\left(T\right)}^{\gamma}.$$

[10]

We have studied the energy and the structure of the OCN (the tree that minimizes Eq.

**10**), first on a general graph and then (in more detail) for the particular case of a discrete lattice (Fig. 4 and*SI Appendix*). Every OCN is a natural river tree, meaning that there exists a height function which determines the tree, and we characterized the natural river trees of a graph in terms of forbidden substructures. In the case of a $d$-dimensional lattice, we determined the energy of the OCN up to a constant factor for every $d\ge 2$ and $\gamma \in \left(\mathrm{0,1}\right)$, extending and improving previous work (8, 12).Fig. 4.

We have also verified that, while it is natural to expect that any OCN on a discrete lattice should take the highly symmetric form of a Peano structure (1, 21), this is in fact not the case: For example, on the square lattice, several suboptimal minimum energy structures are obtained by breaking the Peano structure symmetry (1) (

*SI Appendix*). Much effort has been devoted to reconciling the features of the ground state (known exactly since ref. 10) with those of feasible configurations. The resulting exponents associated with the global minimum did not match either the observational data or the numerical simulations (1). Because every local minimum of the OCN functional is a stationary solution of the general landscape evolution equation, any self-organization of the fluvial landscape corresponds to the dynamical settling of optimal structures into suboptimal niches of their fitness landscape. Thus, feasible optimality (i.e., the search for optima that are accessible to the dynamics given the initial conditions and path obstacles) is the rule in the fluvial landscape, and this might apply to a broad spectrum of interfaces in nature (16). Different network shapes give different values of ${\mathcal{E}}_{\gamma}\left(T\right)$ on a discrete lattice. Such a case imposes sharp conditions on the minimizing structures; in particular, a specific subset of all of the possible trees was selected as the natural rivers restricting the quest of minimum in this set.As stated previously, great interest is placed in OCN on general graphs. In a first application of the theorems presented here, we have addressed the constructability of elevation fields and topographies compatible with the planar imprinting of OCNs (Fig. 4 and

*SI Appendix*). Such results imply the possibility of building replicas of statistically identical matrices for ecological interactions. This is due to the flexibility of random search algorithms in a “greedy” mode (*SI Appendix*), that cannot access the absolute minimum but rather suboptimal configurations that are known to reproduce accurately, differently from the ground state, those of natural fluvial landforms forms (16). Interestingly, an application of limit scaling properties as a test of optimality has been proposed for foodwebs (22) and a direct use of OCN landscapes has been made to explain elevational gradients of biodiversity (17).## Materials and Methods

### Proof of *Theorem 2*.

Let $T$ be an OCN in $G$ with parameter $\gamma $, and suppose (for a contradiction) that $T$ contains a path obstacle. Let $k$ be minimal such that there is a $k$-path obstacle in $T$, so there exist paths ${\eta}_{1},\dots ,{\eta}_{k}$ in $T$, where ${\eta}_{i}$ is ${a}_{i}\to {c}_{i}\left(1\right)\to \cdots \to {c}_{i}\left({\ell}_{i}\right)\to {b}_{i}$, and

$${b}_{1}{a}_{2},{b}_{2}{a}_{3},\dots ,{b}_{k}{a}_{1}\in E\left(G\right)\backslash E\left(T\right).$$

Now, for each $i{\in \mathbb{Z}}_{k}$, define a new tree ${T}_{i}$ by removing the edge ${a}_{i}\to {c}_{i}\left(1\right)$ (or ${a}_{i}\to {b}_{i}$ if ${\ell}_{i}=0$) from $T$, and adding the edge ${a}_{i}\to {b}_{i-1}$ (where ${b}_{0}\u2254{b}_{k}$). If $k=1$, then ${T}_{1}$ is a spanning tree of $G$ with lower energy than $T$ (since ${\ell}_{1}>0$), contradicting our assumption that $T$ is an OCN. Let us therefore assume from now on that $k\ge 2$.

To show that ${T}_{i}$ is a tree, we need to use the minimality of $k$. Indeed, if there is a cycle in ${T}_{i}$ then it must contain the edge ${a}_{i}\to {b}_{i-1}$, and therefore there exists a pathin $T$. But now we can construct a $\left(k-1\right)$-path obstacle in $T$ by replacing the paths ${\eta}_{i-1}$ and ${\eta}_{i}$ by the pathcontradicting the minimality of $k$. Hence, ${T}_{i}$ is a tree, and therefore, since $T$ is an OCN, we have ${\mathcal{E}}_{\gamma}\left(T\right)\le {\mathcal{E}}_{\gamma}\left({T}_{i}\right)$ for each $i{\in \mathbb{Z}}_{k}$.

$${b}_{i-1}\to d\left(1\right)\to \cdots \to d\left(\ell \right)\to {a}_{i},$$

$$\begin{array}{ccc}\hfill {a}_{i-1}\to {c}_{i-1}\left(1\right)\to \hspace{0.17em}\cdots \to {c}_{i-1}\left({\ell}_{i}\right)\to {b}_{i-1}\to d\left(1\right)\to \cdots & \hfill & \hfill \\ \hfill \hspace{0.17em}\cdots \to d\left(\ell \right)\to {a}_{i}\to {c}_{i}\left(1\right)\to \cdots \to {c}_{i}\left({\ell}_{i}\right)\to {b}_{i},& \hfill & \hfill \end{array}$$

We next bound ${\mathcal{E}}_{\gamma}\left({T}_{i}\right)$ from above by introducing a new tree, ${T}_{i}^{\prime}$, by removing the edge ${a}_{i}\to {b}_{i-1}$ from ${T}_{i}$, and adding the edge ${a}_{i}\to {c}_{i-1}\left(1\right)$. Note that this may not be a subgraph of $G$; however, it is a tree [since, as above, there is no path in $T$ from ${c}_{i-1}\left(1\right)$ to ${a}_{i}$, by the minimality of $k$], and therefore its energy ${\mathcal{E}}_{\gamma}\left({T}_{i}^{\prime}\right)$ may still be defined via Eq. For each $i{\in \mathbb{Z}}_{k}$, let us write ${\xi}_{i}$ for the (unique) path in $T$ from ${c}_{i}\left(1\right)$ to the root of $T$. Observe that ${A}_{v}\left({T}_{i}^{\prime}\right)={A}_{v}\left(T\right)-{w}_{i}$ for each $v\in V\left({\xi}_{i}\right)\backslash V\left({\xi}_{i-1}\right)$, where ${w}_{i}\u2254{A}_{{a}_{i}}$, and that ${A}_{v}\left({T}_{i}^{\prime}\right)={A}_{v}\left(T\right)+{w}_{i}$ for each $v\in V\left({\xi}_{i-1}\right)\backslash V\left({\xi}_{i}\right)$, soNow, let us write ${S}_{v}\u2254\left\{i{\in \mathbb{Z}}_{k}:v\in V\left({\xi}_{i}\right)\right\}$ for each vertex $v\in V\left(G\right)$, and define ${V}_{R}\u2254\left\{v\in V\left(G\right):{S}_{v}=R\right\}$ for each set $R{\subset \mathbb{Z}}_{k}$. Note that the vertices in ${V}_{R}$ form a subpath of each ${\xi}_{i}$ such that $i\in R$, and setSince ${\mathcal{E}}_{\gamma}\left(T\right)\le {\mathcal{E}}_{\gamma}\left({T}_{i}^{\prime}\right)$, it follows from Eq. Now, observe that ${f}_{R}$ is a concave function of $x$ for any $\gamma \in \left(\mathrm{0,1}\right)$ and $R{\subset \mathbb{Z}}_{k}$, and note that ${f}_{R}\left(0\right)=0$. It follows that there exists a constant ${\alpha}_{R}\in \mathbb{R}$ such that ${f}_{R}\left(x\right)\le {\alpha}_{R}x$ for every $x\in \mathbb{R}$, and moreover ${f}_{R}\left(x\right)<{\alpha}_{R}x$ whenever $x\ne 0$ and ${V}_{R}\ne \varnothing $. Hence, it follows from Eq. and moreover, the inequality is strict if ${V}_{R}\ne \varnothing $ for any set $R$ included in either sum. Hence, adding the ${\alpha}_{R}$ with $\left\{i-1,i\right\}\subset R$ to both sides of Eq. since ${c}_{i}\left(1\right)\in {V}_{\left\{i\right\}}$, so ${V}_{\left\{i\right\}}\ne \varnothing $. Setting ${\alpha}_{i}\u2254{\sum}_{R:i\in R}{\alpha}_{R}$, it follows that ${\alpha}_{1}>{\alpha}_{2}>\cdots >{\alpha}_{k}>{\alpha}_{1}$, which is the desired contradiction.$\square $

**1**. Moreover, we have ${\mathcal{E}}_{\gamma}\left({T}_{i}\right)\le {\mathcal{E}}_{\gamma}\left({T}_{i}^{\prime}\right)$, since ${A}_{v}\left({T}_{i}\right)\le {A}_{v}\left({T}_{i}^{\prime}\right)$ for every vertex $v\in V\left(G\right)$, and, hence, by the observations above,$${\mathcal{E}}_{\gamma}\left(T\right)\u2a7d{\mathcal{E}}_{\gamma}\left({T}_{i}\right)\u2a7d{\mathcal{E}}_{\gamma}\left({T}_{i}^{\prime}\right).$$

$$\begin{array}{ll}\hfill {\mathcal{E}}_{\gamma}\left({T}_{i}^{\prime}\right)={\hspace{0.17em}\mathcal{E}}_{\gamma}\left(T\right)& +\sum _{v\in V\left({\xi}_{i}\right)\backslash V\left({\xi}_{i-1}\right)}({({A}_{v}\left(T\right)-{w}_{i})}^{\gamma}-{A}_{v}{\left(T\right)}^{\gamma})\hfill \\ \hfill & +\sum _{v\in V\left({\xi}_{i-1}\right)\backslash V\left({\xi}_{i}\right)}({\left({A}_{v}\left(T\right)+{w}_{i}\right)}^{\gamma}-{A}_{v}{\left(T\right)}^{\gamma}).\hfill \end{array}$$

[11]

$${f}_{R}\left(x\right)\u2254\sum _{v\in {V}_{R}}({\left({A}_{v}\left(T\right)+x\right)}^{\gamma}-{A}_{v}{\left(T\right)}^{\gamma}).$$

**11**that$$\sum _{\begin{array}{c}R:i\in R,\\ i-1\notin R\end{array}}{f}_{R}\left(-{w}_{i}\right)+\sum _{\begin{array}{c}R\hspace{0.17em}:\hspace{0.17em}i\notin R,\\ i-1\in R\end{array}}{f}_{R}\left({w}_{i}\right)\ge 0.$$

[12]

**12**, and the fact that ${w}_{i}>0$, that$$\sum _{\begin{array}{c}R:i\in R,\\ i-1\notin R\end{array}}{\alpha}_{R}\le \sum _{\begin{array}{c}R:i\notin R,\\ i-1\in R\end{array}}{\alpha}_{R},$$

[13]

**14**, we obtain$$\sum _{R:i\in R}{\alpha}_{R}<\sum _{R:i-1\in R}{\alpha}_{R},$$

[14]

## Acknowledgments

P.B. and B.B. are partly supported by NSF Grant DMS 1600742. J.B. is partly supported by NSF Grant DMS-1500121 and an Arnold O. Beckman Research Award (University of Illinois Urbana–Champaign Campus Research Board 15006). B.B. and G.C. are partially supported by MULTIPLEX Grant 317532. R. Morris is partially supported by Brazilian National Council for Scientific and Technological Development Grant 303275/2013-8, Fundação Carlos Chagas Filho de Amparo à Pesquisa do Estado do Rio de Janeiro Grant 201.598/2014, and European Research Council (ERC) Starting Grant 680275 MALIG. A.R. and E.B. have been supported by ERC Advanced Grant RINEC 22761. Part of the research in this paper was carried out while J.B. was a Visiting Fellow Commoner at Trinity College, Cambridge.

## Supporting Information

Appendix (PDF)

- Download
- 1.37 MB

## References

1

I Rodríguez-Iturbe, A Rinaldo

*Fractal River Basins: Chance and Self-Organization*(Cambridge Univ Press, Cambridge, UK, 2001).2

B Bollobás

*Graph Theory, An Introductory Course*(Springer, 1st Ed, New York, 1979).3

B Bollobás

*Random Graphs*(Academic, London, 1985).4

I Rodríguez-Iturbe, A Rinaldo, R Rigon, R Bras, E Ijjasz-Vasquez, Fractal structures as least energy patterns: The case of river networks.

*Geophys Res Lett***19**, 889–892 (1992).5

A Rinaldo, et al., Minimum energy and fractal structure of drainage networks.

*Water Resour Res***28**, 2183–2190 (1992).6

J Banavar, F Colaiori, A Flammini, A Maritan, A Rinaldo, Topology of the fittest transportation network.

*Phys Rev Lett***84**, 4745–4748 (2000).7

J Banavar, F Colaiori, A Flammini, A Maritan, A Rinaldo, Scaling, optimality and landscape evolution.

*J Stat Phys***104**, 1–33 (2001).8

LA Briggs, M Krishnamoorthy, Exploring network scaling through variations on optimal channel networks.

*Proc Natl Acad Sci USA***110**, 19295–300 (2013).9

A Rinaldo, I Rodríguez-Iturbe, R Rigon, E Ijjasz-Vasquez, R Bras, Self-organized fractal river networks.

*Phys Rev Lett***70**, 822–825 (1993).10

A Maritan, F Colaiori, A Flammini, M Cieplak, J Banavar, Disorder, river patterns and universality.

*Science***272**, 984–988 (1996).11

A Rinaldo, A Maritan, F Colaiori, A Flammini, R Rigon, Thermodynamics of fractal networks.

*Phys Rev Lett***76**, 3364–3367 (1996).12

F Colaiori, A Flammini, A Maritan, JR Banavar, Analytical and numerical study of optimal channel networks.

*Phys Rev E***55**, 1298–1310 (1997).13

JR Banavar, A Maritan, A Rinaldo, Size and form in efficient transportation networks.

*Nature***399**, 130–132 (1999).14

F Santambrogio, Optimal channel networks, landscape function and branched transport.

*Inter Free Boundaries***9**, 149–169 (2007).15

A Abed-Elmdoust, A Singh, Z Yang, Emergent spectral properties of river network topology: An optimal channel network approach.

*Sci Rep***7**, 11486 (2017).16

A Rinaldo, R Rigon, JR Banavar, A Maritan, I Rodriguez-Iturbe, Evolution and selection of river networks: Statics, dynamics, and complexity.

*Proc Natl Acad Sci USA***111**, 2417–2424 (2014).17

E Bertuzzo, et al., Geomorphic controls on elevational gradients of species richness.

*Proc Natl Acad Sci USA***113**, 1737–1742 (2016).18

A Abed-Elmdoust, M Miri, A Singh, Reorganization of river networks under changing spatiotemporal precipitation patterns: An optimal channel network approach.

*Water Resour Res***52**, 8845–3049 (2016).19

P Bak, C Tang, K Wieselfeld, Self-organized criticality. An explanation of 1/f noise.

*Phys Rev Lett***59**, 381–384 (1987).20

D Dhar, S Majumdar, Abelian sandpile model on Bethe lattice.

*J Phys A: Math Gen***23**, 4333–4350 (1990).21

A Flammini, F Colaiori, Exact analysis of the Peano basin.

*J Phys A Math Gen***29**, 6701–6708 (1996).22

D Garlaschelli, G Caldarelli, L Pietronero, Universal scaling relations in food webs.

*Nature***423**, 165–168 (2003).## Information & Authors

### Information

#### Published in

#### Classifications

#### Copyright

Copyright © 2018 the Author(s). Published by PNAS. This open access article is distributed under Creative Commons Attribution-NonCommercial-NoDerivatives License 4.0 (CC BY-NC-ND).

#### Submission history

**Published online**: June 11, 2018

**Published in issue**: June 26, 2018

#### Keywords

#### Acknowledgments

P.B. and B.B. are partly supported by NSF Grant DMS 1600742. J.B. is partly supported by NSF Grant DMS-1500121 and an Arnold O. Beckman Research Award (University of Illinois Urbana–Champaign Campus Research Board 15006). B.B. and G.C. are partially supported by MULTIPLEX Grant 317532. R. Morris is partially supported by Brazilian National Council for Scientific and Technological Development Grant 303275/2013-8, Fundação Carlos Chagas Filho de Amparo à Pesquisa do Estado do Rio de Janeiro Grant 201.598/2014, and European Research Council (ERC) Starting Grant 680275 MALIG. A.R. and E.B. have been supported by ERC Advanced Grant RINEC 22761. Part of the research in this paper was carried out while J.B. was a Visiting Fellow Commoner at Trinity College, Cambridge.

### Authors

#### Competing Interests

The authors declare no conflict of interest.

## Metrics & Citations

### Metrics

#### Citation statements

#### Altmetrics

### Citations

If you have the appropriate software installed, you can download article citation data to the citation manager of your choice. Simply select your manager software from the list below and click Download.

#### Cited by

Loading...

## View Options

### View options

#### PDF format

Download this article as a PDF file

DOWNLOAD PDF### Get Access

#### Login options

Check if you have access through your login credentials or your institution to get full access on this article.

Personal login Institutional Login#### Recommend to a librarian

Recommend PNAS to a Librarian#### Purchase options

Purchase this article to get full access to it.