River landscapes and optimal channel networks

Balister, Paul; Balogh, József; Bertuzzo, Enrico; Bollobás, Béla; Caldarelli, Guido; Maritan, Amos; Mastrandrea, Rossana; Morris, Robert; Rinaldo, Andrea

doi:10.1073/pnas.1804484115

River landscapes and optimal channel networks

Paul Balister, József Balogh, Enrico Bertuzzo, +5 , Béla Bollobás, Guido Caldarelli, Amos Maritan, Rossana Mastrandrea, Robert Morris, and Andrea Rinaldo [email protected]Authors Info & Affiliations

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

June 11, 2018

115 (26) 6548-6553

https://doi.org/10.1073/pnas.1804484115

PDF/EPUB

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 \geq 2

. Results extend our capabilities in environmental statistical mechanics.

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: (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 (G)

of vertices (or nodes), which correspond to sites in a drainage river basin, and a set

E (G)

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

Δ h_{v} = max_{u \in N (v) \cup {v}} {h_{v} - h_{u}},

where

N (v) = {u : u v \in E (G)}

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 (v) \cup {v}

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} Δ 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).

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 \geq 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 / (d + 1)}

(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

γ \in (0,1)

, we can associate to the (rooted) spanning trees of a graph

G

an energy function of the following form:

E_{γ} (T) ≔ \sum_{v \in V (G)} A_{v} {(T)}^{γ},

[1]

where

A_{v} (T)

is the number of vertices

u \in V (G)

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. 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} (T)

represents the total flow of water out of

v

. The energy function defined in 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:

E_{γ} (T) = \sum_{v \in V (G)} m_{v} g Δ h_{v},

[2]

where

m_{v}

is the mass of water leaving vertex

v

, so

m_{v}

is proportional to

A_{v} (T)

. Furthermore, experimental evidence suggests the following relationship between

A_{v} (T)

and the difference of heights of two adjacent nodes:

Δ h_{v} \propto A_{v} {(T)}^{- 1 / 2}

(1). We therefore obtain:

E_{γ} (T) \propto \sum_{v \in V (G)} A_{v} {(T)}^{1 / 2} .

[3]

Replacing the exponent

1 / 2

with the parameter

γ \in (0,1)

, and ignoring the (unimportant) constant factor, gives Eq. 1. When

γ = 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 (G) \to R

such that if

v \in V (G)

is a nonroot, and

u \in N (v)

, then

v \to u

is an edge of

T

if and only if

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

[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 (v)

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

E_{γ} (T)

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

(η_{1}, \dots, η_{k})

in

T

, with the following property: For each

i {\in Z}_{k}

, the last vertex of

η_{i}

is connected to the first vertex of

η_{i + 1}

by an edge of

G

that is not in

T

.

In other words, if

η_{i}

is the path

a_{i} \to c_{i} (1) \to \dots \to c_{i} (ℓ_{i}) \to b_{i}

in

T

, and

b_{1} a_{2}, b_{2} a_{3}, \dots, b_{k} a_{1} \in E (G) \ E (T),

then we say that

T

contains a

k

-path obstacle (Fig. 1). We emphasize that the paths

η_{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

ℓ_{1} > 0

.

Fig. 1.

k-path obstacles. (A) A three-path obstacle. For $i \in {1,2,3}$ , $η_{i}$ , is the directed path $a_{i} \to c_{i} (1) \to \dots \to c_{i} (ℓ_{i}) \to b_{i}$ in $T$ , where $ℓ_{1}, ℓ_{2}, ℓ_{3}$ are arbitrary nonnegative integers. Dashed red lines represent the edges $b_{1} a_{2}$ , $b_{2} a_{3}$ and $b_{3} a_{1}$ belonging to $G$ but not to $T$ . (B) An example of a two-path obstacle ( $a_{1} \to c_{1} \to b_{1}$ and $a_{2} \to c_{2} \to b_{2}$ ), where $G$ (dashed red lines) is a regular 2D lattice and $T$ (black arrows) is a spanning tree.

It is not difficult to show that a natural river tree cannot contain a

k

-path obstacle for any

k \in 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 \geq 1

, a natural river tree

T

cannot contain a

k

-path obstacle. Indeed, let

f : V (G) \to R

be the height function associated to

T

, and suppose that there exist paths

η_{1}, \dots, η_{k}

in

T

, where

η_{i}

is the path

a_{i} \to c_{i} (1) \to \dots \to c_{i} (ℓ_{i}) \to b_{i}

in

T

, and

b_{1} a_{2}, b_{2} a_{3}, \dots, b_{k} a_{1} \in E (G) \ E (T)

. Since

a_{i} \to c_{i} (1)

(or

a_{i} \to b_{i}

if

ℓ_{i} = 0

) is in

T

, and

a_{i} b_{i - 1} \in E (G)

, it follows from the definition that

f (b_{i}) \leq f (c_{i} (1)) < f (b_{i - 1})

for each

i {\in Z}_{k}

, where

b_{0} ≔ b_{k}

. But then,

f (b_{1}) > f (b_{2}) > \dots > f (b_{k}) > f (b_{1}),

[5]

and this contradiction proves the claim.

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 (G)

and

\begin{matrix} E (D) = {x \to y : either x \to y is an edge of T, or \\ \exists z such that x z \in E (G) \ E (T) and (z \to y) \in E (T)} . \end{matrix}

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 \geq 1

, since

T

is a tree. Let us denote the vertices of the cycle as follows:

\begin{matrix} c_{1} (1) \to \dots \to c_{1} (ℓ_{1}) \to b_{1} \to c_{2} (1) \to \dots \to c_{2} (ℓ_{2}) \to \\ b_{2} \to \dots \to b_{k - 1} \to c_{k} (1) \to \dots \to c_{k} (ℓ_{k}) \to b_{k} \to c_{1} (1), \end{matrix}

where

b_{1}, \dots, b_{k}

are the vertices whose out-neighbor in the cycle is not an edge of

T

, and

ℓ_{1}, \dots, ℓ_{k} \geq 0

. Thus, there exist vertices

a_{1}, \dots, a_{k}

such that, for each

i {\in Z}_{k}

, we have

b_{i - 1} a_{i} \in E (G) \ E (T)

(where

b_{0} ≔ b_{k}

), and

a_{i} \to c_{i} (1)

(or

a_{i} \to b_{i}

if

ℓ_{i} = 0

) is an edge of

T

. This gives a

k

-path obstacle, contradicting our assumption.

□

It follows from the claim that the directed graph

D

can be extended to a linear ordering

<

on

V (G)

(i.e., a linear ordering with

x > y

for every edge

x \to y

of

D

). Choose a height function

f

such that

f (x) < f (y)

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 (G)

be a nonroot, and note that there is a unique

u \in N (v)

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 (u) < f (v)

. Now let

u \neq w \in N (v)

, and observe that

w \to u

is an edge of

D

(since

w v \in E (G) \ E (T)

and

v \to u

is an edge of

T

), so

f (u) < f (w)

. Hence,

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

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 Theorem 1.

□

Recall that an OCN (with parameter

γ

) in a graph

G

is a spanning tree of

G

that minimizes the energy

E_{γ} (T)

. 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 N

and

γ \in (0,1)

. If

T

is an OCN in

G

with parameter

γ

, 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.

In other words, if

G

is a finite graph with a single root, and

T

is an OCN in

G

with parameter

γ \in (0,1)

, then

T

contains no path obstacle, by Theorem 2, and hence, by Theorem 1, is a natural river spanning tree of

G

.

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

{[n]}^{d}

, with root

(1, \dots, 1)

, and edges between vertices at

ℓ_{\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 N

and

γ \in (0,1)

, define

T (n, d, γ) ≔ min {_{γ} (T) : T is a spanning tree of {[n]}^{d}},

the energy of an OCN in the

d

-dimensional grid. The following theorem determines

T (n, d, γ)

up to a constant factor (depending on

d

and

γ

) for all values of the parameters. It extends results of Colaiori et al. (12), who proved the lower bounds for

d = 2

and

γ \neq 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 (n)

in the case

d = 2

and

γ = 1 / 2

.

Theorem 4.

For every integer

d \geq 2

, we have

T (n, d, γ) = \{\begin{matrix} Θ (n^{d}), & i f 0 < γ < 1 - 1 / d, \\ Θ (n^{d} \log n), & i f γ = 1 - 1 / d, \\ Θ (n^{1 + d γ}), & i f 1 - 1 / d < γ < 1, \end{matrix}

[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

{[n / 2]}^{d}

, oriented so that

2^{d} - 1

of them have their root near the center of

{[n]}^{d}

, and the last one has its root at

(1, \dots, 1)

, Fig. 2. Join the vertex

(n / 2, \dots, n / 2)

in this last subcube to each of the roots of the other subcubes, which are each of the form

(n / 2, \dots, n / 2) + x

for some

x \in {0,1}^{d}

, so as to form a spanning tree

T_{n}

in

{[n]}^{d}

.

Fig. 2.

The upper bound construction in the proof of *Theorem 4*. The tree is formed (in two dimensions) by four copies of the construction for ${[n / 2]}^{2}$ , placed so that squares 1, 2, and 3 have their root adjacent to the center of ${[n]}^{2}$ , and the last one has its root at (1,1) (the black dot in the left bottom corner).

We claim that the subtree

T_{n}

of

{[n]}^{d}

constructed in this way contains the diagonal edge

(x + 1, \dots, x + 1) \to (x, \dots, x)

for each

1 \leq 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

(n / 2 + 1, \dots, n / 2 + 1) \to (n / 2, \dots, n / 2)

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

{[n / 2]}^{d}

flows down the diagonal of the final copy, and hence we have the following recursive formula:

\begin{matrix} E_{γ} (T_{n}) & = & 2^{d} \cdot E_{γ} (T_{n / 2}) + \\ + \sum_{x = 1}^{n / 2} {((A_{(x, \dots, x)} (T_{n / 2}) + n^{d} - {(n / 2)}^{d}}^{γ} + \\ - A_{(x, \dots, x)} {(T_{n / 2})}^{γ}) . \end{matrix}

[7]

Since

{(y + z)}^{γ} \leq y^{γ} + z^{γ}

, and noting that

(n / 2) {(n^{d} - {(n / 2)}^{d}}^{γ} \leq n^{1 + d γ}

, it follows that

E_{γ} (T_{n}) \leq 2^{d} \cdot E_{γ} (T_{n / 2}) + n^{1 + d γ},

and, hence, noting that

E_{γ} (T_{1}) = 1

, and recalling that

n = 2^{k}

, we obtain

T (n, d, γ) \leq E_{γ} (T_{n}) \leq n^{1 + d γ} + \dots + 2^{ℓ d} {(\frac{n}{2^{ℓ}}}^{1 + d γ} + \dots + n^{d} .

[8]

This inequality immediately implies the claimed upper bounds. Indeed, if

1 + d γ < d

then the right-hand side of Eq. 8 is an increasing geometric series, so

T (n, d, γ) ⩽ \frac{n^{d}}{1 - 2^{1 + d γ - d}} = O (n^{d}) .

If

1 + d γ > d

then it is a decreasing geometric series, so

T (n, d, γ) ⩽ \frac{n^{1 + d γ}}{1 - 2^{d - 1 - d γ}} = O (n^{1 + d γ}) .

Finally, if

1 + d γ = d

then all terms are equal, so

T (n, d, γ) ⩽ n^{d} \log_{2} n + 1 = O (n^{d} \log n),

as required.

We next turn to the lower bounds. Observe first that the lower bound

E_{γ} (T) = \sum_{v \in V ({[n]}^{d})} A_{v} {(T)}^{γ} ⩾ \sum_{v \in V ({[n]}^{d})} 1 = n^{d},

holds trivially for every spanning tree

T

of

{[n]}^{d}

, and every

γ \in (0,1)

, since we have

A_{v} (T) \geq 1

for every vertex

v

. Hence,

T (n, d, γ) \geq n^{d}

holds always.

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

{[n]}_{○}^{d}

denote a copy of the

d

-dimensional grid in which every vertex on the boundary is a root and define

\begin{array}{l} F (n, d, γ) ≔ min {\{E}_{γ} (F) : F is a rooted \\ spanning forest of {[n]}_{○}^{d}}, \end{array}

where each component of

F

is assumed to have a single root (which lies on the boundary of

{[n]}^{d}

), and the energy

E_{γ} (F)

is defined to be the sum of the energies of the components. Note that

F (n, d, γ) \leq T (n, d, γ)

, since a spanning tree in

{[n]}^{d}

contains a rooted spanning forest of

{[n]}_{○}^{d}

. We claim that

F (3 n, d, γ) \geq 3^{d} \cdot F (n, d, γ) + \frac{γ}{3^{d (1 - γ)}} \cdot n^{1 + d γ} .

[9]

To prove Eq. 9, we partition

{[3 n]}^{d}

into

3^{d}

subcubes of size

{[n]}^{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 (n, d, γ)

to the energy. Additionally, the middle subcube has to send its mass from its boundary to the boundary of the large cube

{[3 n]}^{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} (T) \leq {(3 n)}^{d}

and

γ < 1

, we have

A_{v} {(T)}^{γ} - {(A_{v} (T) - x}^{γ} \geq γ \cdot A_{v} {(T)}^{γ - 1} \cdot x \geq γ \cdot {(3 n)}^{d γ - d} \cdot x,

for any

x > 0

. Therefore, moving the mass of the center subcube to the boundary increases the energy by at least

n^{d} \cdot n \cdot γ \cdot {(3 n)}^{d γ - d} = \frac{γ}{3^{d (1 - γ)}} \cdot n^{1 + d γ},

and, hence, we obtain Eq. 9 (Fig. 3).

Fig. 3.

The proof of the lower bounds in *Theorem 4*. Black dots represent the root on the boundary of ${[3 n]}^{d}$ for each component of the spanning forest. The red square highlights the central subcube; the water from this subcube must travel to a root on the boundary, and to do so it must first cross its own boundary and then pass through at least $n$ additional vertices.

Setting

C = γ \cdot 3^{- (d + 1)}

, and noting that

F (1, d, γ) = 1

, we obtain

\begin{matrix} T (n, d, γ) \geq F (n, d, γ) \geq C \cdot n^{1 + d γ} + \dots \\ \dots + C \cdot 3^{ℓ d} {(\frac{n}{3^{ℓ}}}^{1 + d γ} + \dots + n^{d}, \end{matrix}

where in the penultimate term

ℓ = k - 2

. If

1 + d γ < d

then this gives only a minor improvement over the trivial lower bound

n^{d}

, but for

1 + d γ > d

the improvement is more substantial. Indeed, it implies that

T (n, d, γ) ⩾ \frac{γ + o (1)}{3^{d + 1}} \cdot \frac{n^{1 + d γ}}{1 - 3^{d - 1 - d γ}} = Ω n^{1 + d γ},

as

n \to \infty

. Finally, if

1 + d γ = d

, then we obtain

T (n, d, γ) ⩾ \frac{γ}{3^{d + 1}} \cdot n^{d} \log_{3} n - 2 = Ω n^{d} \log n,

as required.

□

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 (1)

, and hence to determine

T (n, d, γ)

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

E_{γ} (T) \propto \sum_{v \in V (G)} A_{v} {(T)}^{γ} .

[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 \geq 2

and

γ \in (0,1)

, extending and improving previous work (8, 12).

Fig. 4.

OCN and landscapes. (A) An OCN derived within a 64 $\times$ 64 regular 2D lattice by means of a simulated annealing approach (*SI Appendix*) assuming $γ = 1 / 2$ . The algorithm ensures that the state is a local, dynamically accessible minimum, not necessarily the absolute minimum of the function $E_{γ} (T)$ . (B) The corresponding landscape (i.e., the elevation $h_{v}$ ) computed assuming $Δ h_{v} \propto A_{v} {(T)}^{- 1 / 2}$ . We have verified numerically that the obtained OCN is a natural river spanning tree.

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

E_{γ} (T)

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

γ

, 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

η_{1}, \dots, η_{k}

in

T

, where

η_{i}

is

a_{i} \to c_{i} (1) \to \dots \to c_{i} (ℓ_{i}) \to b_{i}

, and

b_{1} a_{2}, b_{2} a_{3}, \dots, b_{k} a_{1} \in E (G) \ E (T) .

Now, for each

i {\in Z}_{k}

, define a new tree

T_{i}

by removing the edge

a_{i} \to c_{i} (1)

(or

a_{i} \to b_{i}

if

ℓ_{i} = 0

) from

T

, and adding the edge

a_{i} \to b_{i - 1}

(where

b_{0} ≔ b_{k}

). If

k = 1

, then

T_{1}

is a spanning tree of

G

with lower energy than

T

(since

ℓ_{1} > 0

), contradicting our assumption that

T

is an OCN. Let us therefore assume from now on that

k \geq 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 path

b_{i - 1} \to d (1) \to \dots \to d (ℓ) \to a_{i},

in

T

. But now we can construct a

(k - 1)

-path obstacle in

T

by replacing the paths

η_{i - 1}

and

η_{i}

by the path

\begin{matrix} a_{i - 1} \to c_{i - 1} (1) \to \dots \to c_{i - 1} (ℓ_{i}) \to b_{i - 1} \to d (1) \to \dots \\ \dots \to d (ℓ) \to a_{i} \to c_{i} (1) \to \dots \to c_{i} (ℓ_{i}) \to b_{i}, \end{matrix}

contradicting the minimality of

k

. Hence,

T_{i}

is a tree, and therefore, since

T

is an OCN, we have

E_{γ} (T) \leq E_{γ} (T_{i})

for each

i {\in Z}_{k}

.

We next bound

E_{γ} (T_{i})

from above by introducing a new tree,

T_{i}^{'}

, by removing the edge

a_{i} \to b_{i - 1}

from

T_{i}

, and adding the edge

a_{i} \to c_{i - 1} (1)

. 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} (1)

to

a_{i}

, by the minimality of

k

], and therefore its energy

E_{γ} (T_{i}^{'})

may still be defined via Eq. 1. Moreover, we have

E_{γ} (T_{i}) \leq E_{γ} (T_{i}^{'})

, since

A_{v} (T_{i}) \leq A_{v} (T_{i}^{'})

for every vertex

v \in V (G)

, and, hence, by the observations above,

E_{γ} (T) ⩽ E_{γ} (T_{i}) ⩽ E_{γ} (T_{i}^{'}) .

For each

i {\in Z}_{k}

, let us write

ξ_{i}

for the (unique) path in

T

from

c_{i} (1)

to the root of

T

. Observe that

A_{v} (T_{i}^{'}) = A_{v} (T) - w_{i}

for each

v \in V (ξ_{i}) \ V (ξ_{i - 1})

, where

w_{i} ≔ A_{a_{i}}

, and that

A_{v} (T_{i}^{'}) = A_{v} (T) + w_{i}

for each

v \in V (ξ_{i - 1}) \ V (ξ_{i})

, so

\begin{align} E_{γ} (T_{i}^{'}) = E_{γ} (T) & + \sum_{v \in V (ξ_{i}) \ V (ξ_{i - 1})} ({(A_{v} (T) - w_{i})}^{γ} - A_{v} {(T)}^{γ}) \\ + \sum_{v \in V (ξ_{i - 1}) \ V (ξ_{i})} ({(A_{v} (T) + w_{i})}^{γ} - A_{v} {(T)}^{γ}) . \end{align}

[11]

Now, let us write

S_{v} ≔ {i {\in Z}_{k} : v \in V (ξ_{i})}

for each vertex

v \in V (G)

, and define

V_{R} ≔ {v \in V (G) : S_{v} = R}

for each set

R {\subset Z}_{k}

. Note that the vertices in

V_{R}

form a subpath of each

ξ_{i}

such that

i \in R

, and set

f_{R} (x) ≔ \sum_{v \in V_{R}} ({(A_{v} (T) + x)}^{γ} - A_{v} {(T)}^{γ}) .

Since

E_{γ} (T) \leq E_{γ} (T_{i}^{'})

, it follows from Eq. 11 that

\sum_{\begin{array}{c} R : i \in R, \\ i - 1 \notin R \end{array}} f_{R} (- w_{i}) + \sum_{\begin{array}{c} R : i \notin R, \\ i - 1 \in R \end{array}} f_{R} (w_{i}) \geq 0 .

[12]

Now, observe that

f_{R}

is a concave function of

x

for any

γ \in (0,1)

and

R {\subset Z}_{k}

, and note that

f_{R} (0) = 0

. It follows that there exists a constant

α_{R} \in R

such that

f_{R} (x) \leq α_{R} x

for every

x \in R

, and moreover

f_{R} (x) < α_{R} x

whenever

x \neq 0

and

V_{R} \neq \emptyset

. Hence, it follows from Eq. 12, and the fact that

w_{i} > 0

, that

\sum_{\begin{array}{c} R : i \in R, \\ i - 1 \notin R \end{array}} α_{R} \leq \sum_{\begin{array}{c} R : i \notin R, \\ i - 1 \in R \end{array}} α_{R},

[13]

and moreover, the inequality is strict if

V_{R} \neq \emptyset

for any set

R

included in either sum. Hence, adding the

α_{R}

with

{i - 1, i} \subset R

to both sides of Eq. 14, we obtain

\sum_{R : i \in R} α_{R} < \sum_{R : i - 1 \in R} α_{R},

[14]

since

c_{i} (1) \in V_{{i}}

, so

V_{{i}} \neq \emptyset

. Setting

α_{i} ≔ \sum_{R : i \in R} α_{R}

, it follows that

α_{1} > α_{2} > \dots > α_{k} > α_{1}

, which is the desired contradiction.

□

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).

Google Scholar

2

B Bollobás Graph Theory, An Introductory Course (Springer, 1st Ed, New York, 1979).

Google Scholar

3

B Bollobás Random Graphs (Academic, London, 1985).

Google Scholar

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).

Crossref

Google Scholar

5

A Rinaldo, et al., Minimum energy and fractal structure of drainage networks. Water Resour Res 28, 2183–2190 (1992).

Crossref

Google Scholar

6

J Banavar, F Colaiori, A Flammini, A Maritan, A Rinaldo, Topology of the fittest transportation network. Phys Rev Lett 84, 4745–4748 (2000).

Crossref

PubMed

Google Scholar

7

J Banavar, F Colaiori, A Flammini, A Maritan, A Rinaldo, Scaling, optimality and landscape evolution. J Stat Phys 104, 1–33 (2001).

Crossref

Google Scholar

8

LA Briggs, M Krishnamoorthy, Exploring network scaling through variations on optimal channel networks. Proc Natl Acad Sci USA 110, 19295–300 (2013).

Crossref

PubMed

Google Scholar

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).

Crossref

PubMed

Google Scholar

10

A Maritan, F Colaiori, A Flammini, M Cieplak, J Banavar, Disorder, river patterns and universality. Science 272, 984–988 (1996).

Crossref

PubMed

Google Scholar

11

A Rinaldo, A Maritan, F Colaiori, A Flammini, R Rigon, Thermodynamics of fractal networks. Phys Rev Lett 76, 3364–3367 (1996).

Crossref

PubMed

Google Scholar

12

F Colaiori, A Flammini, A Maritan, JR Banavar, Analytical and numerical study of optimal channel networks. Phys Rev E 55, 1298–1310 (1997).

Crossref

Google Scholar

13

JR Banavar, A Maritan, A Rinaldo, Size and form in efficient transportation networks. Nature 399, 130–132 (1999).

Crossref

PubMed

Google Scholar

14

F Santambrogio, Optimal channel networks, landscape function and branched transport. Inter Free Boundaries 9, 149–169 (2007).

Crossref

Google Scholar

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).

Crossref

PubMed

Google Scholar

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).

Crossref

PubMed

Google Scholar

17

E Bertuzzo, et al., Geomorphic controls on elevational gradients of species richness. Proc Natl Acad Sci USA 113, 1737–1742 (2016).

Crossref

PubMed

Google Scholar

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).

Crossref

Google Scholar

19

P Bak, C Tang, K Wieselfeld, Self-organized criticality. An explanation of 1/f noise. Phys Rev Lett 59, 381–384 (1987).

Crossref

PubMed

Google Scholar

20

D Dhar, S Majumdar, Abelian sandpile model on Bethe lattice. J Phys A: Math Gen 23, 4333–4350 (1990).

Crossref

Google Scholar

21

A Flammini, F Colaiori, Exact analysis of the Peano basin. J Phys A Math Gen 29, 6701–6708 (1996).

Crossref

Google Scholar

22

D Garlaschelli, G Caldarelli, L Pietronero, Universal scaling relations in food webs. Nature 423, 165–168 (2003).

Crossref

PubMed

Google Scholar

Information & Authors

Information

Published in

Proceedings of the National Academy of Sciences

Vol. 115 | No. 26
June 26, 2018

PubMed: 29891709

Classifications

Copyright

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

Affiliations

Paul Balister

Department of Mathematical Sciences, University of Memphis, Memphis, TN 38152;

View all articles by this author

József Balogh

Department of Mathematics, University of Illinois, Urbana, IL 61801;

View all articles by this author

Enrico Bertuzzo

Dipartimento di Scienze Ambientali e Informatica, Università Cà Foscari Venezia, 30123 Venice, Italy;

View all articles by this author

Béla Bollobás

Department of Mathematical Sciences, University of Memphis, Memphis, TN 38152;

Trinity College, Cambridge CB2 1TQ, United Kingdom;

London Institute for Mathematical Sciences, London W1K 2XF, United Kingdom;

View all articles by this author

Guido Caldarelli

Networks Unit, IMT (Istituzioni, Mercati, Tecnologie) School for Advanced Studies Lucca, 55100 Lucca, Italy;

Institute of Complex Systems, Consiglio Nazionale delle Ricerche, 00185 Rome, Italy;

European Centre for Living Technology, Università Cà Foscari Venezia, 30124 Venice, Italy;

View all articles by this author

Amos Maritan

Dipartimento di Fisica e Astronomia, University of Padua, 35122 Padua, Italy;

Istituto Nazionale di Fisica Nucleare, University of Padua, 35122 Padua, Italy;

View all articles by this author

Rossana Mastrandrea

Networks Unit, IMT (Istituzioni, Mercati, Tecnologie) School for Advanced Studies Lucca, 55100 Lucca, Italy;

View all articles by this author

Robert Morris

Instituto de Matemática Pura e Aplicada, Jardim Botânico, Rio de Janeiro, RJ 22460-320, Brazil;

View all articles by this author

Andrea Rinaldo¹ [email protected]

Laboratory of Ecohydrology, École Polytechnique Fédérale de Lausanne, Lausanne CH-1015, Switzerland;

Dipartimento di Ingegneria Civile, Edile e Ambientale, University of Padua, 35122 Padua, Italy

View all articles by this author

Notes

1

To whom correspondence should be addressed. Email: [email protected].

Author contributions: P.B., J.B., B.B., and R. Morris obtained the mathematical proofs; and P.B., J.B., E.B., B.B., G.C., A.M., R. Mastrandrea, R. Morris, and A.R. determined the problem, discussed the results, and wrote the paper.

Reviewers: A.A.-E., University of Texas at Austin; A.S., ETH Zurich; and R.V., Microsoft Research Redmond.

Competing Interests

The authors declare no conflict of interest.

Metrics & Citations

Metrics

Note: The article usage is presented with a three- to four-day delay and will update daily once available. Due to ths delay, usage data will not appear immediately following publication. Citation information is sourced from Crossref Cited-by service.

Altmetrics

Citations

Cite this article

P. Balister,
J. Balogh,
E. Bertuzzo,
B. Bollobás,
G. Caldarelli,
A. Maritan,
R. Mastrandrea,
R. Morris,
& A. Rinaldo,

River landscapes and optimal channel networks, Proc. Natl. Acad. Sci. U.S.A. 115 (26) 6548-6553, https://doi.org/10.1073/pnas.1804484115 (2018).

Copied!

Copying failed.

Export the article citation data by selecting a format from the list below and clicking Export.

Cited by

View Options

View options

PDF format

Download this article as a PDF file

DOWNLOAD PDF

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

Save for later

Purchase options

Purchase this article to access the full text.

Single Article Purchase

River landscapes and optimal channel networks

Proceedings of the National Academy of Sciences

Vol. 115
No. 26
pp. 6513-E6098

Restore content access

Restore content access for purchases made as a guest

Featured Topics

Articles By Topic

Featured Topics

Articles By Topic

Featured Topic

Articles By Topic

Significance

Abstract

Sign up for PNAS alerts.

Definition.

Which Networks Are River Trees?

Definition.

Theorem 1.

Proof.

Claim.

Proof.

Theorem 2.

Theorem 3.

Energy of OCNs in d-Dimensional Grids.

Theorem 4.

Proof.

Discussion

Materials and Methods

Proof of Theorem 2.

Acknowledgments

Supporting Information

References

Information

Published in

Classifications

Copyright

Submission history

Keywords

Acknowledgments

Authors

Affiliations

Notes

Competing Interests

Metrics

Altmetrics

Citations

Cite this article

Cited by

View options

PDF format

Login options

Recommend to a librarian

Purchase options

Restore content access

Share

Share article link

Share on social media

Further reading in this issue

Biological tissue-inspired tunable photonic fluid

Origin of slow spontaneous resting-state neuronal fluctuations in brain networks

YES1 amplification is a mechanism of acquired resistance to EGFR inhibitors identified by transposon mutagenesis and clinical genomics

Bodily maps of emotions

Neural network retuning and neural predictors of learning success associated with cello training

Oxytocin-enforced norm compliance reduces xenophobic outgroup rejection

Sign up for thePNAS Highlights newsletter

Energy of OCNs in $d$ -Dimensional Grids.