††thanks: This work is supported by NSF award DMS-1412557.††thanks: This work
is supported by NSF award EAR-1723033.††thanks: This work is supported by NSF
awards DMS-1839336 and EAR-1811909.

# Critical Tokunaga model for river networks

Yevgeniy Kovchegov<EMAIL_ADDRESS>Department of Mathematics,
Oregon State University, Corvallis, OR, 97331-4605, USA  Ilya Zaliapin
<EMAIL_ADDRESS>Department of Mathematics and Statistics, University of Nevada,
Reno, NV, 89557-0084, USA  Efi Foufoula-Georgiou<EMAIL_ADDRESS>Departments of
Civil and Environmental Engineering and Earth System Science,
University of California, Irvine, CA 92697.

###### Abstract

The hierarchical organization and self-similarity in river basins have been
topics of extensive research in hydrology and geomorphology starting with the
pioneering work of Horton in 1945. Despite significant theoretical and applied
advances however, the mathematical origin of and relation among Horton laws
for different stream attributes remain unsettled. Here we capitalize on a
recently developed theory of random self-similar trees to introduce a one-
parametric family of self-similar critical Tokunaga trees that elucidates the
origin of Horton laws, Hack’s laws, basin fractal dimension, power-law
distributions of link attributes, and power-law relations between distinct
attributes. The proposed family includes the celebrated Shreve’s random
topology model and extends to trees that approximate the observed river
networks with realistic exponents. The results offer tools to increase our
understanding of landscape organization under different hydroclimatic
forcings, and to extend scaling relationships useful for hydrologic prediction
to resolutions higher that those observed.

##  I Introduction

In a pioneering study “of streams and their drainage basins”, Robert E. Horton
Horton45  introduced the concept of river stream order and formulated two
fundamental laws of the composition of stream-drainage nets. The Law of Stream
Numbers postulates a geometric decay of the numbers $N_{K}$ of streams of
increasing order $K$, with the exponent $R_{B}$. The Law of Stream Lengths
postulates a geometric growth of the average length $L_{K}$ of streams of
increasing order $K$, with the exponent $R_{L}$. During the 20th century,
geometric dependence on the stream order (conventionally referred to as Horton
law) has been documented for multiple stream attributes, including upstream
area, magnitude (number of upstream headwater channels, also called sources),
the total channel length, the longest stream length, link slope, mean annual
discharge, energy expenditure, etc. RIR01 . Despite their elemental role in
describing the key regularities in river stream networks (such as fractal
dimension, Hack’s law, etc.), Horton’s laws remain an empirical finding and
their origin and apparent ubiquity remain unsettled Kirchner93 .

The first attempt at a rigorous explanation of Horton laws was made by Ronald
L. Shreve Shreve66  in the 1960s, who examined a “topologically random
population of channel networks”, where all topologically distinct networks
with given number of first-order streams are equally likely. This model is
equivalent to the celebrated critical binary Galton-Watson branching process
with a given population size BWW00 ; Pitman . Shreve’s calculations imply that
in this model the Horton law of stream numbers holds with $R_{B}=4$. Although
not attempted by Shreve, it can be shown BWW00  that the law of stream lengths
also holds here with $R_{L}=2$ under the assumption of constant or identically
distributed edge lengths. Albeit insightful and mathematically tractable, the
random topology model deviates from observations, which became apparent with
the development of improved methods for extracting river networks DBE94 ;
Pec95 . This called for developing alternative modeling approaches for river
networks.

It has proven challenging to find a model that would be mathematically
tractable and flexible enough to reproduce the Horton exponents and other
scaling laws observed in river basins. One end of the modeling spectrum is
occupied by conceptual approaches, such as Peano fractal basin (RIR01 , their
Sect. 2.4; Horton45 , their Fig. 25) or Scheidegger’s lattice model Sch67 ;
TNT88 . These models provide an invaluable insight into the origin of the
observed scalings; they however lack realistic dendritic patterns and values
of scaling exponents. On the other end are simulation approaches that can
generate visually appealing networks that closely fit selected exponents, but
can be analytically opaque. The Optimal Channel Network (OCN) model RRR+92 ;
RRR+93 ; RRRIB93 ; MRR+96 ; RRBMR14 ; BBB+18  is a well-recognized simulation
technique. Following the energy expenditure minimization principle, the model
creates random drainage basins on a planar lattice (or more general graphs).
We refer to RIR01  for a comprehensive discussion of these and other models.

Despite the progress achieved by the modeling efforts of the $20^{th}$
century, the following essential questions remain unanswered: What are
sufficient conditions for Horton laws? What are the values of the Horton
exponents for river basins? How are the Horton exponents for different stream
attributes related to each other and to other basin parameters? There is a
consensus that Horton’s laws are connected to the self-similar structure of a
basin RIR01 ; Pec95 ; GW89 ; GTD07 , which is generally understood as
invariance of basin’s statistical structure under changing the scale of
analysis (zooming in or out) Turcotte_book . However, a consensus is still
lacking about a suitable mathematical definition of tree self-similarity.
Three alternative definitions have been proposed: the Toeplitz property of the
Tokunaga coefficients Pec95 ; NTG97 ; the invariance of a tree distribution
with respect to the Horton pruning (cutting source streams) BWW00 ; and
statistical self-similarity of basin attributes GW89 ; PG99 . The related
unsettled questions are: How are the alternative definitions related? and Is
self-similarity (any version) sufficient for selected Horton laws? These
questions extend to other areas beyond hydrogeomorphology where Horton laws
and related scalings have been reported, including computer science F+79 ;
DP06 , statistical seismology BP04 ; HTR08 ; Y13 ; ZBZ13 , vascular analysis
Kassab00 , brain studies Cetal06 , ecology CLF07 , and biology TPN98 .

We answer the above questions within a self-consistent mathematical theory of
random self-similar trees recently developed by the authors KZ20survey . The
theory develops the concept of tree self-similarity that unifies the
alternative existing definitions, rigorously explains the appearance and
parameterization of Horton laws, and offers a novel approach to modeling a
variety of dendritic systems. The goal of this paper is to adapt and extend
the theory to the studies of river networks. Most notably, we show that two
fundamental and practically appealing properties – coordination and Horton
prune invariance – result in trees that respect a wealth of scaling laws
observed in landscape dissection. Furthermore, we propose a new one-parametric
family of critical Tokunaga trees, which reproduces multiple Horton laws and
related scalings reported in river network studies, with realistic values of
the respective parameters. The critical Tokunaga family yields rigorous
relations among scaling exponents that have been empirically documented in
multiple earlier studies, and serves as a useful analytic and modeling tool
for further analysis of river network structure and dynamics. Our results also
offer a computationally efficient algorithm of generating self-similar trees
with arbitrary parameters (Horton exponents, fractal dimensions, etc.), which
facilitates ensemble simulations.

We represent a stream network that drains a single basin (watershed,
catchment) as a rooted binary tree. The basin outlet (point furthest
downstream) corresponds to the tree root, sources (points furthest upstream)
to leaves, junctions (points where two streams meet) to internal vertices, and
links (stream segments between two successive nodes) to edges. This graph-
theoretic nomenclature provides a link to the probability and combinatorics
literature on the topic. We assume that all examined trees belong to the space
${\mathcal{L}}$ of finite binary rooted planted trees with positive edge
lengths KZ20survey . Recall that a tree is called planted if the tree root has
degree $1$ (the most downstream link goes to an ocean or another large water
body instead of merging with another stream). The space ${\mathcal{L}}$
includes the empty tree $\phi$ comprised of the root vertex. We also consider
the space ${\mathcal{T}}$ of combinatorial projections of trees from
${\mathcal{L}}$, that is trees with the same combinatorial structure but no
edge lengths.

##  II Results

###  II.1 Review of Horton laws and implied scaling relationships

The Horton-Strahler orders for river streams have been discussed extensively
in the literature and nicely reviewed and summarized in RIR01 . In this
section we introduce the orders through the viewpoint of Horton pruning – this
streamlines our exposition and prepares the reader for material that follows.
For completeness, we also briefly recall Horton laws and their implications.

Consider the map ${\mathcal{R}}:{\mathcal{L}}\to{\mathcal{L}}$ that removes
the source links from a tree ${\mathcal{T}}$. The Horton-Strahler order of a
tree ${\mathcal{T}}$ Horton45 ; Strahler ; Melton59  is the minimal number of
Horton prunings that completely eliminates it:

${\sf ord}(T)=\min\left\\{k\geq 0\,:\,{\mathcal{R}}^{k}(T)=\phi\right\\}.$ (1)

The Horton pruning and Horton-Strahler orders are illustrated in Fig. 1 (see
Appendix IV.1 for details and an alternative computational definition).

Figure 1:  Horton pruning of the stream network of Beaver creek, Floyd County,
KY. Streams of orders $K=2,\ldots,6$ are shown by different colors (see legend
on the right). Streams of order $1$ (source streams) are not shown for visual
convenience. (A) The first Horton pruning, after eliminating streams of order
K=1. (B)-(E) Second to fifth consecutive Horton prunings. The sixth pruning
completely eliminates the network. The channel extraction is done using
RiverTools software (http://rivix.com).

Horton’s Law of Stream Numbers Horton45  postulates a geometric decay of the
stream counts $N_{K}$ of increasing order $K$ with exponent $R_{B}\geq 2$:

$N_{K}~{}\propto R_{B}^{-K}\quad\text{ or }\quad\frac{N_{K}}{N_{K+1}}=R_{B}.$
(2)

The notation $x\propto y$ stands for $x={\rm Const.}\times y$. The lower bound
on $R_{B}$ follows from the definition of Horton-Strahler orders, since it
takes at least two streams of order $K$ to create a single stream of order
$K+1$. The value of $R_{B}$ in large river basins is usually close to $4.5$
RIR01 ; Shreve66 ; GW89 ; Turcotte_book ; PG99 ; Strahler ; Leopold ;
Tarboton96 ; GW98 ; ZZF13 ; Mesa18 . Figure 2A (cyan circles) shows the Horton
law for stream numbers in the stream network of Beaver Creek of Fig. 1; here
$R_{B}\approx 4.6$.

Figure 2:  Critical Tokunaga fit to the Horton laws in the stream network of
Beaver creek. The stream network is shown in Fig. 1. Symbols correspond to the
values of the observed attributes. Lines and dots show the fit by the critical
Tokunaga model with $c=2.3$. (A) Stream numbers $N_{K}$. The model fit is
given by (13); it has asymptotic slope $-\log_{10}(2c)\approx-0.66$. (B)
Average stream magnitude $M_{K}$ (cyan circles) and average number of links
$S_{K}$ (magenta squares). The fit for $M_{K}$ is given by (14); it has
asymptotic slope $\log_{10}(2c)\approx 0.66$. The fit for $S_{K}$ is given by
$c^{K-1}$; it has theoretical slope $\log_{10}c\approx 0.36$. (C) Average
total contributing area $A_{K}$ (cyan circles) and average total upstream
channel length $L_{K}^{\rm tot}$ (magenta squares). The fitting lines,
according to a combination of equations (14) and (50), have theoretical slope
$\log_{10}(2c)\approx 0.66$. (D) Average stream length $L_{K}$ (magenta
squares) and average length $\Lambda_{K}$ of the longest stream to the divide
(cyan circles). The fitting lines, according to equations (21) and (23), have
theoretical slope $\log_{10}(c)\approx 0.36$.

Horton’s Law of Stream Lengths Horton45  postulates a geometric growth of the
average length $L_{K}$ of streams of increasing order $K$ with exponent
$R_{L}$:

$L_{K}~{}\propto R_{L}^{K}\quad\text{ or }\quad{L_{K+1}\over L_{K}}=R_{L}.$
(3)

The value of $R_{L}$ in river networks is around $2.5$ RIR01 . Figure 2D
(magenta squares) shows the Horton law for stream lengths in the Beaver Creek
network of Fig. 1; here $R_{L}\approx 2.3$.

Similarly to the laws (2),(3) discussed above, a geometric scaling of any
average stream attribute $Z_{K}$ with order $K$, is also called Horton law,
and the respective exponent is called Horton exponent Pec95 ; NTG97 . Horton
laws are documented for multiple physical and combinatorial attributes,
including upstream area, magnitude (number of upstream sources), total
upstream channel length, length of the longest channel to the divide, etc.
RIR01 . Figure 2 shows Horton laws for seven stream attributes of the Beaver
Creek that is shown in Fig. 1. We use a convention that the Horton exponent is
greater than unity, which always can be achieved by selecting the sign of the
exponent K in the Horton law (e.g., Eqs. (2),(3)).

Horton laws imply power-law frequencies of link attributes and power-law
relations between the average values of distinct attributes. Suppose that
stream attributes $Z$ and $Y$ satisfy Horton’s law with Horton exponents
$R_{Z}$ and $R_{Y}$, respectively. Using each of the laws to express the
channel order $K$ and equating these expressions, we find

$Z_{K}~{}\propto~{}Y_{K}^{\alpha},\quad\text{ with }\quad\alpha=\frac{\log
R_{Z}}{\log R_{Y}}.$ (4)

Equation (4) is a punctured (by discrete orders) version of a power-law
relation $Z\propto Y^{\alpha}$ that is abound among hydrologic quantities. A
well-studied example is Hack’s law that relates the length $L$ of the longest
stream in a basin to the basin area $A$ via $L\propto A^{\bf h}$ with ${\bf
h}\approx 0.6$ Hack57 ; Rigon+96 . Equation (4) suggests that the parameter
${\bf h}$ is expressed via the exponents $R_{L}$ and $R_{A}$ of the Horton
laws for length $L$ and area $A$ as:

${\bf h}=\frac{\log R_{L}}{\log R_{A}}.$ (5)

Next, consider the value $Z_{(i)}$ of an attribute $Z$ calculated at link $i$
in a large basin. Assuming Horton laws for $N_{K}$ and $S_{K}$ with exponents
$R_{B}$ and $R_{S}$, respectively, and considering a limit of an infinitely
large basin we approximate the distribution of link attributes as (see
Appendix IV.3)

$\\#\\{i:Z_{(i)}\geq z\\}~{}\propto~{}z^{-\beta},\quad\beta=\frac{\log
R_{B}-\log R_{S}}{\log R_{Z}}.$ (6)

Such power laws are documented for the upstream contributing area, length of
the longest channel to the divide, water discharge, or energy expenditure. For
example, analyses of river basins (e.g., MRR+96 ; TBR88 ; RI92 ) extracted
from digital elevation models (DEM’s) suggest

$\\#\\{i:A_{(i)}\geq a\\}\propto a^{-\beta_{A}},\quad\beta_{A}\approx 0.43$
(7)

and

$\\#\\{i:\Lambda_{(i)}\geq l\\}\propto
l^{-\beta_{\Lambda}},\quad\beta_{\Lambda}\approx 0.8,$ (8)

where $A_{(i)}$ is the area upstream of link $i$ and $\Lambda_{(i)}$ is the
distance from link $i$ to the furthest source (or, equivalently, to the basin
divide) measured along the channel network.

Horton laws (e.g. Eqs. (2),(3)) and the implied scaling relations (Eqs.
(4),(6)) provide key observational constraints in modeling river stream
networks RIR01 ; RRR+93 ; Mesa18 ; TBR88 . Our work explains the appearance of
Horton laws in terms of tree self-similarity and offers a parametric toolbox
for the analysis and modeling of river networks and other branching structures
that exhibit such scaling relations.

###  II.2 Tree self-similarity and Tokunaga sequence

We introduce the concept of self-similarity for random trees that encompasses
the existing definitions and satisfies practical intuitive expectations. The
proposed definition applies to a distribution of trees from a suitable space
such as ${\mathcal{T}}$ or ${\mathcal{L}}$ and combines two fundamental
properties – coordination and Horton prune-invariance.

Coordination means that the (random) structure of a river basin is determined
by its order. For example, a basin with outlet of order three and a sub-basin
of order three within a basin of order nine have, statistically, the same
structure. This assumption is at the heart of analyses based on the Horton-
Strahler orders; it has been imposed, explicitly or implicitly, in the
mainstream studies of river networks Horton45 ; RIR01 ; Shreve66 ; DBE94 ;
Pec95 ; PG99 ; Tarboton96 . A distribution that satisfies the coordination
property is called coordinated. We refer to KZ20survey  for a measure-
theoretic definition of coordination.

Horton prune-invariance formalizes the expectation that the scaling laws of
hydrology are (by and large) independent of the data resolution. The Horton
pruning ${\mathcal{R}}$ is a natural model for the change of resolution in a
river network. Indeed, better observations lead to detecting smaller streams,
which increases the basin order. Pruning a basin by order is roughly
equivalent to decreasing the resolution of stream detection. The Horton prune-
invariance requires that the statistical structure of trees remains the same
after zooming in or out.

######  Definition 1 (Self-similar tree).

A coordinated distribution $\mu$ on the space $\mathcal{T}$ of combinatorial
trees is called self-similar if it is invariant with respect to Horton pruning
BWW00 ; KZ16 :

$\mu\left(\mathcal{R}^{-1}(T)|T\neq\phi\right)=\mu(T)\quad\text{for any}\quad
T\in\mathcal{T}.$ (9)

Recall that $\phi$ denotes an empty tree. Equation (9) states that, for any
non-empty tree $T$, the total probability assigned by $\mu$ to all trees that
result in $T$ after pruning – these trees are collectively denoted by
$R^{-1}(T)$ – is the same as the probability of $T$. Informally, consider a
forest of trees generated by measure $\mu$, where each tree $T$ occurs
multiple times according to its probability $\mu(T)$. The forest is self-
similar if after pruning each tree by ${\mathcal{R}}$ we obtain the same
forest. This definition can be extended to trees with edge lengths from space
${\mathcal{L}}$; see Def. 9 in KZ20survey . In that case, we allow the edge
lengths to scale after pruning by a multiplicative scaling constant $\zeta>0$.
We use a conventional abuse of terminology by saying that a tree $T$ is self-
similar meaning that $T$ is a random tree drawn from a self-similar
distribution $\mu$.

The measure-theoretic Definition 1 might be not appealing for practical
analyses that oftentimes involve only a handful of finite basins. A bridge
from this definition to easily computed stream statistics is provided by
Theorem 6 (see Appendix IV.4). It shows that the expected value $T_{i,j}$ of
the number $n_{i,j}$ of streams of order $i$ that merge with a randomly
selected stream of order $j$ in a self-similar tree is a function of the
difference $j-i$ (and not individual values of $i$ and $j$). Accordingly, each
self-similar measure $\mu$ corresponds to a unique non-negative sequence of
Tokunaga coefficients

$T_{i,i+k}=T_{k}={\sf E}_{\mu}[n_{i,i+k}]\quad\text{for all}\qquad i,k>0,$
(10)

where ${\sf E}_{\mu}$ denotes the mathematical expectation with respect to
$\mu$. Moreover, we show (see Cor. 3 in Appendix IV.4) that the Toeplitz
property $T_{i,i+k}=T_{k}$ and Horton prune-invariance of Eq. (9) are
equivalent for coordinated measures (i.e., both hold or do not hold at the
same time). Hence, the Tokunaga coefficients $T_{k}$ provide a fundamental
parameterization of a self-similar tree and constitute the main tool of
respective applied analyses.

Our Definition 1 of tree self-similarity unifies the alternative definitions
used in the literature. Burd et al. BWW00  define self-similarity in Galton-
Watson trees as the Horton prune-invariance; this is a special case of our
definition since the Galton-Watson trees are coordinated KZ20survey . Peckham
Pec95  and Newman et al. NTG97  define self-similarity as the Toeplitz
property for Tokunaga coefficients; this is equivalent to our definition in
coordinated trees (Corollary 3). The coordination assumption is further
justified in KZ20survey  by showing that the Toeplitz property alone, without
coordination, allows for a multitude of obscure measures that are hardly
useful in practice. Gupta and Waymire GW89  and Peckham and Gupta PG99
suggested a concept of statistical self-similarity that requires a random
stream attribute $Z$ to have distribution that scales with order. It can be
shown (Sect. 7 in KZ20survey ) that (i) statistical self-similarity for some
attributes (e.g., for any discrete attribute) may only hold asymptotically,
and (ii) multiple attributes, including stream length, magnitude, and total
basin length, are statistically self-similar in a limit of infinitely large
basin that is self-similar according to our Definition 1.

###  II.3 Horton laws for stream numbers, magnitudes in self-similar trees

We now capitalize on the concept of tree self-similarity introduced above to
establish a key emergent property of self-similar trees – Horton laws for
stream numbers and magnitudes, conveniently parameterized by the Tokunaga
sequence.

Consider the mean number

${\mathcal{N}}_{i}[K]={\sf E}_{\mu}[N_{i}(T)\,|\,{\sf ord}(T)=K]$

of streams of order $i$ in a basin of order $K$ and the mean magnitude (number
of upstream sources) $M_{i}$ of a stream of order $i$. Consider also the
generating function $T(z)=\sum\limits_{k=1}^{\infty}T_{k}z^{k}$ of the
Tokunaga coefficients and define

$\hat{t}(z)=-1+2z+T(z).$ (11)

######  Theorem 1 (Horton law for stream numbers, magnitudes in self-similar
tree KZ20survey ; KZ16 ).

Consider a self-similar tree $T$ with Tokunaga sequence $T_{k}$ and suppose
that

$\limsup\limits_{k\to\infty}T_{k}^{1/k}<\infty.$ (12)

Then, the stream numbers ${\mathcal{N}}_{i}[K]$ and the magnitudes $M_{i}$
obey Horton laws:

$\lim_{K\to\infty}\left({\mathcal{N}}_{1}[K]\,R_{B}^{-K}\right)=M~{}<\infty,$
(13) $\lim_{i\to\infty}\left(M_{i}\,R_{M}^{-i}\right)=M~{}<\infty.$ (14)

The Horton exponents are given by

$R_{B}=R_{M}=1/w_{0},$ (15)

where $w_{0}$ is the only real root of the function $\hat{t}(z)$ in the
interval $(0,1/2]$ and $M$ is a positive real constant given by

$M=-\frac{1}{w_{0}}\lim_{z\to w_{0}}\frac{z(z-w_{0})}{\hat{t}(z)}.$ (16)

The proof is given in Appendix IV.5.

Theorem 1 shows that the Horton laws for mean stream numbers
${\mathcal{N}}_{i}[K]$ and magnitudes $M_{i}$ hold in almost any self-similar
tree, more specifically – in a tree with an arbitrary Tokunaga sequence
$T_{k}$. The only restriction of Eq. (12) prohibits super-exponential growth
of $T_{k}$, such as $T_{k}=k!$ or $T_{k}=k^{k}$. The theorem establishes a
strong form of the Horton law (Eqs. (13),(14)), which implies a weaker version
that is often reported in applied literature:

$\displaystyle\lim\limits_{K\to\infty}\frac{{\mathcal{N}}_{i}[K]}{{\mathcal{N}}_{i+1}[K]}$
$\displaystyle=$ $\displaystyle R_{B}~{}~{}~{}\text{for any }i,\text{ and }$
$\displaystyle\lim\limits_{i\to\infty}\frac{M_{i+1}}{M_{i}}$ $\displaystyle=$
$\displaystyle R_{M}.$ (17)

Theorem 1 emphasizes the existence of a multitude of self-similar measures
with the same Horton exponent. Assume we fix $R_{B}$ and hence the root
$w_{0}$ of $\hat{t}(z)$ according to Eq. (15). Equation (11) readily asserts
that there is an infinite number of Tokunaga sequences that correspond to an
arbitrary $w_{0}$ within $(0,1/2]$. For example, if $R_{B}=4$, then
$w_{0}=1/4$ and one needs $T(z)=1/2$. This can be achieved by selecting any of
$T_{k}=\\{2,0,\ldots\\}$, $\,\\{1,4,0,\ldots\\}$, $\,\\{0,8,0,\ldots\\}$,
$\,\\{2^{k-1}\\}$, etc., where “$\ldots$” denotes trailing zeros.

We observe that the Horton law of stream numbers in Theorem 1 (Eqs. (13)) is
an asymptotic statement, different from the ideal Horton law for stream
numbers (2) which is commonly used in the literature. This is not a
mathematical peculiarity – the ideal Horton law is merely an approximation to
the actual behavior of stream counts. Its approximate nature is not related to
the finite size of the observed basins – the ideal Horton law rarely holds in
theoretical trees of arbitrarily large size. Formally, we show in Appendix
IV.6 that the ideal Horton law for stream numbers in a self-similar tree holds
if and only if $T_{k}=0$ for $k>1$. Realistically, Horton laws are asymptotic
statements of different strength. For example, the strongest form of Horton
law for stream numbers is that of Eq. (13), which implies a weaker version of
Eq. (II.3). Accordingly, the power relations among different stream attributes
(4) and power-law frequencies of link attributes (6) that we have derived from
the ideal Horton law of Eq. (2) remain heuristic. A formal analysis based on
actual Horton laws (like those in Eqs. (13) and (14)), which will be presented
elsewhere, confirms the results of Eqs. (4) and (6) and reveals additional
solutions with oscillatory tail behavior.

The Horton laws for other stream attributes may or may not hold depending on
additional assumptions about $T_{k}$ and other details of basin organization.
A comprehensive treatment is possible using the generating function approach
outlined in Appendix IV.2. Most importantly, further analysis often requires
specifying a concrete self-similar distribution, not only its Tokunaga
sequence $T_{k}$. Below we examine a particularly useful family of
distributions.

###  II.4 Random Attachment Model (RAM) of self-similar trees

According to Theorem 6 (Appendix IV.4), every self-similar measure corresponds
to a unique Tokunaga sequence $T_{k}$. At the same time, a multitude of self-
similar measures can be constructed for a given Tokunaga sequence. Here we
introduce a particularly symmetric random tree (tree distribution) for a given
Tokunaga sequence and establish its key properties. We use Poisson attachment
construction within exponential segments; this ensures that the link lengths
have exponential distribution and the attachment of streams of lower orders to
a given stream of a larger order is done in uniform random fashion. We refer
to this construction as Random Attachment Model (RAM).

The RAM specifies a tree distribution on ${\mathcal{L}}$ by a non-negative
Tokunaga sequence $T_{k}$, the order distribution $\pi_{K}={\sf P}({\sf
ord}(T)=K)$, and the distribution of stream lengths. The model assumes that
the lengths of streams of order $j$ are independent exponential random
variables with rate $\lambda_{j}$. Hence, the model is specified by three non-
negative sequences:

$\\{T_{k}\\},~{}~{}\\{\lambda_{j}\\},~{}\text{ and }~{}\\{\pi_{K}\\}.$

Importantly, the probabilities $\pi_{K}$ and rates $\lambda_{j}$ do not affect
the combinatorial structure of a tree of a given order, which is completely
specified by the Tokunaga coefficients $T_{k}$.

A random tree is constructed in a hierarchical fashion, starting from the
stream of the highest order and adding side-tributaries of consecutively
smaller orders. The tree order $K$ is selected according to the distribution
$\pi_{K}$. At the first step we generate the main stream that will have order
$K$ in the final tree; its length is an exponential random variable with rate
$\lambda_{K}$. At each of the remaining $K-1$ steps, we add streams of lower
orders to the existing tree by a Poisson attachment procedure. The streams
added at step $m$ will have order $i(m)=K-m+1$ in the final tree. The lengths
of the newly added streams are independent exponential random variables with
rate $\lambda_{i(m)}$. The new streams are added in two steps. First, we
consider the existing tree as a one-dimensional metric space (union of link
segments) and generate a collection of points on this space according to a
homogeneous Poisson process. The process intensity depends on the order of a
link within the final tree. Specifically, within every link added at step
$K-j+1$ (that will have order $j$ in the final tree) the Poisson intensity is
$\lambda_{j}T_{j-i(m)}$. A single stream is then attached to each Poisson
point. Second, we add two new streams to each source stream of the current
tree (except the sources just added during this step). The first part of this
procedure (Poisson attachment) ensures that the tree has Tokunaga coefficients
$T_{k}$, and the second part (adding stream pairs) increases the tree order by
one at each step.

Figure 3: Examples of critical Tokunaga trees. The trees are generated by the
critical Tokunaga process with $c=2.3$ and Horton-Strahler order $K=5$. The
line width is proportional to the contributing area. The figure accurately
represents the tree combinatorial structure; the edge lengths are scaled for a
better planar embedding.

The trees generated by RAM can be equivalently represented as trajectories of
a continuous-time multitype Hierarchical Branching Process, with time evolving
from the root upstream and member types corresponding to the stream orders.
This approach, explored by the authors in KZ20survey , yields the joint
distribution of the orders $K_{1}<K_{2}$ of subtrees that share a common root
of order $K$:

${\sf P}\left(K_{1}=j,K_{2}=m|K\right)=\begin{cases}S^{-1}_{K}&\\!\\!\text{if
}j=m=K-1,\\\ T_{K-j}S^{-1}_{K}&\\!\\!\text{if }j<m=K.\end{cases}$ (18)

We now use this result to propose a computationally efficient recursive
construction of RAM trees. A tree of order $1$ is a stream of exponential
length with rate $\gamma$. To create a tree of order $K>1$ we first generate a
link of exponential length with rate $\lambda_{K}S_{K}$, where

$S_{K}=1+\sum_{i=1}^{K-1}T_{i}.$ (19)

To this link we attach two conditionally (conditioned on the order $K$)
independent trees whose orders are drawn from Eq.(18). Each of these trees is
generated according to the same recursive procedure. This algorithm generates
trees with up to $10^{6}$ edges within seconds, providing a flexible
computational framework for ensemble simulations based on independent
statistical realizations of a tree with fixed parameters. Examples of RAM
stream networks are shown in Fig. 3.

Another useful result of the branching process theory establishes the
necessary and sufficient conditions for a RAM tree to be self-similar
according to Definition 1. These conditions explicitly parameterize the
probabilities $\pi_{K}$ and stream length rates $\lambda_{j}$ for an arbitrary
Tokunaga sequence $T_{k}$. This emphasizes the richness of self-similar
family.

######  Theorem 2 (Self-similar RAM; Thm. 11 in KZ20survey ).

Suppose $T$ is a random tree generated by the RAM with parameters
$\\{T_{k}\\}$, $\\{\lambda_{j}\\}$, and $\\{\pi_{K}\\}$. Then $T$ is a
coordinated tree. Tree $T$ is self-similar with scaling constant $\zeta>0$
(see Definition 1 and its discussion) if and only if

$\displaystyle\pi_{K}$ $\displaystyle=$ $\displaystyle
p(1-p)^{K-1}~{}~{}(K\geq 1)~{}~{}~{}\text{ and }$ $\displaystyle\lambda_{j}$
$\displaystyle=$ $\displaystyle\gamma\,\zeta^{1-j}~{}~{}~{}(j\geq 1)$ (20)

for some parameters $p\in(0,1)$, $\gamma>0$ and $\zeta>0$ (and any Tokunaga
sequence $T_{k}$).

######  Corollary 1 (Horton law for the stream lengths).

Consider a self-similar RAM tree with parameters given by Eq. (2). Then the
average length $L_{j}$ of a stream of order $j$ satisfies

$L_{j}R_{L}^{-j}=\frac{1}{\zeta\gamma}<\infty\quad\text{ with }\quad
R_{L}=\zeta.$ (21)

The proof is given in Appendix IV.8. We notice that the Horton law holds here
in an exact form, without a limit in order $j$.

We show below that the two well-known properties of self-similar river basins
– fractal dimension and Horton law for the longest stream length – formally
hold in a self-similar RAM tree.

######  Theorem 3 (Fractal dimension of a self-similar tree).

Consider a self-similar RAM tree with Tokunaga sequence $T_{k}$ and other
parameters given by Eq. (2). Let $w_{0}$ be the only real root of the function
$\hat{t}(z)$ (Eq. (11)) in the interval $(0,1/2]$. Then, the fractal dimension
of the tree in the limit of infinite order and after a suitable length
rescaling, is given by

${\bf d}=\max\\{1,{\bf d}_{0}\\},\quad{\bf d}_{0}=-\frac{\log
w_{0}}{\log\zeta}=\frac{\log R_{B}}{\log R_{L}}.$ (22)

The proof is given in Appendix IV.9.

Equation (22) coincides with the expression first obtained by La Barbera and
Rosso LR89  using a heuristic assumption of a basin with an ideal Horton law
of stream numbers. Figure 4A shows a map of ${\bf d}$ as a function of the
Horton exponents $R_{B}$ and $R_{L}$.

Figure 4:  Fractal dimension ${\bf d}$ (panel A) of equation (22) and Hack’s
exponent ${\bf h}={\bf d}^{-1}$ (panel B) of a self-similar RAM tree in the
limit of infinite size as a function of the Horton exponents $R_{B}$ and
$R_{L}$. Selected levels of ${\bf d}$ and ${\bf h}$ are shown by marked black
lines. Green thick lines correspond to the critical Tokunaga process
(Definition 2, Eq. (24)) for which $R_{B}=2R_{L}$. Blue dots depict the pairs
$(R_{B},R_{L})$ estimated in nine real river basins by TBR88 , see also
(RIR01, , Table 2.1).

######  Theorem 4 (Horton law for the length of the longest stream).

Consider a self-similar RAM basin with parameters given by Eq. (2). Let
$\Lambda_{k}$ denote the average length of the longest stream in a basin of
order $k$. Then

$\lim\limits_{k\to\infty}\Lambda_{k}R_{\Lambda}^{-k}={\rm
Const.}<\infty\quad\text{ with }\quad R_{\Lambda}=\zeta.$ (23)

The proof is given in Appendix IV.10.

###  II.5 Critical Tokunaga tree and emergent scaling relations

Observations on river networks have supported a basic constraint that the link
length distribution is independent of the position of the link within a basin
RIR01 ; TBR89 . This motivates one to describe a family of trees that respect
this property. Surprisingly, this leads to a one-parameter family of critical
Tokunaga trees that satisfy multiple additional symmetries and include the
celebrated Shreve model as a special case.

The length of a link of order $K$ in the RAM model is an exponential random
variable with rate $\lambda_{K}S_{K}$, which is a direct consequence of using
Poisson attachment along exponentially distributed streams. The order-
independent link length implies $\lambda_{K}S_{K}={\rm Const.}$ Using the
general form of $\lambda_{K}$ in a self-similar RAM tree (2), one can select
$\gamma$ such that $S_{K}=\zeta^{K-1}$. This corresponds to the unique form of
the Tokunaga sequence in a self-similar RAM tree with identically distributed
link lengths:

$T_{k}=(c-1)\,c^{k-1}~{}~{}~{}(k\geq 1),\qquad c=\zeta\geq 1.$ (24)

######  Definition 2 (Critical Tokunaga tree).

A self-similar RAM with $p={1\over 2}$ and Tokunaga sequence of Eq. (24) is
called critical Tokunaga tree or critical Tokunaga model. By Theorem 2, for
this model we have $\pi_{K}=2^{-K}$ and $\lambda_{j}=\gamma\,\zeta^{1-j}$.

Figure 3 shows three examples of critical Tokunaga trees of order $K=5$ with
parameter $c=2.3$, which gives a close approximation to the observed river
networks (see Table 1). Equation (24) is a special case of the two-parameter
sequence $T_{k}=ac^{k-1}$ introduced by Eliji Tokunaga Tok78  to approximate
river basin branching; this sequence has been examined in detail in Pec95 ;
Turcotte_book ; NTG97 ; ZZF13 ; KZ16 . The one-parameter sequence of Eq. (24)
appears in the random self-similar network (RSN) model of Veitzer and Gupta
VG00 , which uses a purely combinatorial algorithm of recursive local
replacement of the network generators to construct random trees. A theoretical
underpinning of this constraint is revealed via the prism of branching process
analysis. Kovchegov and Zaliapin KZ18  have shown that a random tree $T$
generated by the critical Tokunaga model is critical and time-invariant in
both combinatorial and metric forms KZ20survey . In particular, the condition
$p=1/2$ is necessary and sufficient for criticality. Moreover, the Geometric
Branching Process (that generates the combinatorial part of a RAM tree) is
time invariant if and only if it corresponds to the critical Tokunaga model.
Recall that criticality means that a branching process has unit average
population size after an arbitrary but fixed time advancement (in both
discrete and continuous versions). Time-invariance means that the frequency of
orders of subtrees that survive after a given time advancement coincides with
the initial order distribution.

It is natural to assume that the local contributing area of a link (area that
contributes to the link directly, and not via its descendant joint) is a
function of the link length. This allows us to examine the average total
contributing area $A_{i}$ of a stream of order $i$. In particular, the order-
independent link lengths imply order-independent local areas. The following
result establishes Hack’s law in a critical Tokunaga tree.

######  Theorem 5 (Hack’s law in a critical Tokunaga tree).

Consider a critical Tokunaga tree (Definition 2). Then the average lengths
$\Lambda_{i}$ of the longest stream and the average total contributing area
$A_{i}$ of a basin are related as

$\displaystyle\Lambda_{i}$ $\displaystyle\sim$ $\displaystyle{\rm
Const.}\times\big{(}A_{i}\big{)}^{\\!{\bf h}},\text{ where }$
$\displaystyle{\bf h}$ $\displaystyle=$ $\displaystyle{\bf
d}^{-1}\\!\\!=-{\log\zeta\over\log w_{0}}=\frac{\log R_{L}}{\log R_{B}}.$ (25)

The proof is given in Appendix IV.8.

The Hack’s law of Eq. (5) also holds in more general self-similar RAM trees
(which may not be critical Tokunaga) as is shown in Appendix IV.9. Figure 4B
shows a map of ${\bf h}$ as a function of the Horton exponents $R_{B}$ and
$R_{L}$. The critical Tokunaga case corresponds to $R_{B}=2R_{L}=2c$ and hence
${\bf d}=\log_{c}(2c)$ and ${\bf h}=\log_{2c}c$; this case is depicted by
green line in Fig. 4.

Combining our results, we obtain the following summary of the Horton exponents
in a critical Tokunaga tree.

######  Corollary 2 (Horton exponents in a critical Tokunaga tree).

The Horton exponents in the critical Tokunaga model are given by

$2c=R_{B}=R_{M}=R_{A}>R_{S}=R_{L}=R_{\Lambda}=c.$ (26)

The proof is given in Appendix IV.8.

Corollary 2 reveals that the essential Horton exponents in critical Tokunaga
trees only assume two distinct values ($c$ and $2c$). The inequality
$R_{S}<R_{B}$, which is a part of Eq. (26), has been conjectured by Peckham
Pec95  for trees with a well-defined Tokunaga sequence.

Notably, the critical Galton-Watson process with exponential edge lengths
Pitman , which is equivalent to Shreve’s random topology model after
conditioning on the basin magnitude, is a special case of critical Tokunaga
model with $c=2$ (KZ20survey, , Thm. 15). In other words, the critical
Tokunaga model offers a natural extension of the critical binary Galton-Watson
process to a similarly versatile family of processes with a wide range of
Horton exponents, fractal dimensions, Hack’s exponents, and other parameters.
As such, this model may be useful for multiple fields beyond hydrology.

Results of Chunikhina EVC18 ; EVCthesis  imply that the critical Tokunaga
model with $c=2$ maximizes the entropy rate among the trees that satisfy the
Horton law of stream numbers, and that the critical Tokunaga model with a
fixed $c$ maximizes the entropy rate among the trees that satisfy the Horton
law for stream numbers with $R_{B}=2c$.

Some additional scaling properties of the critical Tokunaga trees are
collected in Appendix IV.7.

###  II.6 Critical Tokunaga model closely fits observations

The critical Tokunaga model provides a very close fit to the data and scaling
relations reported in river studies over the past decades. Table 1 summarizes
the values of the key scaling exponents in the critical Tokunaga model and
compares them with exponents in river network observations and the well-
established OCN model RRR+93 . The table uses the results of Corollary 2, Eq.
(26) (Horton exponents), Theorem 3, Eq. (22) (fractal dimension ${\bf d}$),
Theorem 5, Eq. (5) (Hack’s exponent ${\bf h}$), and Eq. (6) (scaling exponents
$\beta_{\Lambda}$ and $\beta_{A}$).

The critical Tokunaga model fit to the observed data is further illustrated in
the Beaver Creek basin of Fig. 1. Figure 2 shows seven Horton laws fit by the
critical Tokunaga model with $c=2.3$. Specifically, we consider the following
stream attributes averaged over streams of order $K=1,\ldots,6$: the stream
number $N_{K}$ (panel A), the average magnitude $M_{K}$ and the average number
$S_{K}$ of links in a stream (panel B), the average total contributing area
$A_{K}$ and the average total upstream channel length $L^{tot}_{K}$ (panel C),
and the average stream length $L_{K}$ and the average length $\Lambda_{K}$ of
the longest channel to the divide (panel D). The fitting lines correspond to
the critical Tokunaga model predictions, which impressively agree with
observations of all examined stream attributes (see figure caption).

##  III Discussion

A solid body of observational, modeling, and theoretical studies connects
Horton laws, power-law distributions of and power-law relations among river
stream attributes to the self-similar structure of stream networks RIR01 ;
BWW00 ; Pec95 ; GW89 ; GTD07 ; Turcotte_book ; PG99 ; Tarboton96 ; GW98 ;
Mesa18 ; VG00 ; DR00 ; GCO96 ; Cetal98 ; PT00 . We suggest a rigorous
treatment of the appearance and parameterization of Horton laws in river
networks using a recently formulated theory of random self-similar trees
KZ20survey . The proposed framework unifies the existing results and
contributes to explaining the ubiquity of Horton laws in dendritic systems of
arbitrary origin.

The main technical contribution of our work is a rigorous treatment of the
appearance of Horton laws in self-similar trees (Eqs.
(13),(14),(II.3),(21),(23)). We show that the two fundamental properties –
coordination and Horton-prune invariance – necessarily lead to the Horton laws
for stream numbers and magnitudes (Theorem 1). Additional mild assumptions,
like those in the Random Attachment Model (RAM), yield the Horton laws for
multiple other attributes (Theorem 4; Corollaries 1,2), which in turn imply
basin fractal dimension (Theorem 3), Hack’s law (Theorem 5), and other power-
law scaling relations (Eqs. (4),(6)). Our results can be easily extended to
other stream attributes such as stream slope, width, depth and velocity, which
are known to be proportional to a power of the upstream magnitude RIR01 ; GW89
. The developed framework may also facilitate analysis of the width function
LF07  or scaling of hydrologic fluxes GTD07 ; GW90  in self-similar basins.
Such analyses can be done either analytically, or using ensemble simulation
that is facilitated by a fast simulation algorithm for RAM trees.

The self-similarity is defined here (Definition 1) as invariance of a
coordinated tree distribution with respect to the operation of Horton pruning,
which is in accord with the empirical and modeling evidence of the past
decades BWW00 ; Pec95 ; NTG97 ; Melton59 ; KZ16 ; VG00 . This approach unifies
three seemingly distinct definitions of self-similarity that existed in the
literature BWW00 ; Pec95 ; GW89 ; NTG97 . Importantly, each self-similar tree
distribution corresponds to a unique Tokunaga sequence $T_{k}$ that quantifies
merging of branches of distinct orders (Theorem 6, Corollary 3). This provides
a fundamental connection between an abstract measure-theoretic prune-
invariance property and the tangible Tokunaga coefficients that can be
statistically estimated in a single tree.

The family of self-similar distributions (Definition 1) rigorously reproduces
the key geomorphic scalings discovered and reconfirmed during the past 80
years for river basins and summarized by RIR01 ; MRR+96 ; GTD07 ;
Turcotte_book ; DR00 , with a close fit to the examined exponents (Table 1,
Fig. 2). Interestingly, this fit is achieved within a one-parameter family of
critical Tokunaga trees (Definition 2, Eq. (24)). Although trees that satisfy
Eq. (24) (and commonly referred to as Tokunaga trees) have been known for a
long time Pec95 ; Tok78 ; VG00 , only very recently a rigorous understanding
has been gained of the theoretical importance of this constraint within the
general framework of branching processes KZ20survey ; KZ16 ; KZ18 . In
addition, neither the order distribution nor link lengths distribution
(specified by $\pi_{K}$ and $\lambda_{j}$ of Eq. (2)) have been examined in
Tokunaga trees. This justifies our reference to the critical Tokunaga tree as
a novel model. The critical Tokunaga model provides a natural parametric
extension of the critical binary Galton-Watson branching process (and includes
it as a special case with $c=2$), which proved to be an indispensable model in
many areas and remains at the forefront of theoretical and applied research
nearly 150 years after its discovery GW . This hints at deep and not fully
understood symmetries in the structure of river networks. The theory of random
self-similar trees explains the mathematical origin of these symmetries and
offers tools for future exploration.

The presented results might advance applied statistical analysis of river
stream attributes, via mapping all quantities of interest to a single master
parameter $c$ of Eq. (24). Statistical estimation of this parameter can be
designed more effectively than that for a range of distinct yet possibly
related quantities (e.g. Horton exponents). This in turn facilitates global
mapping of river network features and studying possible effects of
hydroclimatic variables on landscape dissection. Corollary 2 shows that
multiple Horton laws examined in this work hold with only two distinct Horton
exponents: $R_{B}=R_{M}=R_{A}=2c$ and $R_{L}=R_{S}=R_{\Lambda}=c$. This
substantial reduction of observed quantities is well supported by data (Table
1, Fig. 2) and might inform a range of modeling and theoretical efforts.

The critical Tokunaga model presents an ultimately symmetric class of trees
characterized by coordination, Horton prune-invariance, criticality, time-
invariance, and identically distributed link lengths (and hence local
contributing areas). Despite these multiple constraints, this class is
surprisingly rich, extending from perfect binary trees ($c=1$) to the famous
Shreve’s random topology model ($c=2$) to the structures reminiscent of the
observed river networks ($c\approx 2.3$) and beyond. While offering a
convenient theoretical and modeling paradigm, the critical Tokunaga model is
merely a subclass of a much broader family of self-similar trees that might
better accommodate for various problem-specific data features. For instance,
Fig. 4 suggests that the observed stream networks tend to cluster around the
critical Tokunaga line $R_{B}=2R_{L}$ in the ($R_{B}$,$R_{L}$) space. An
applied study can use the self-similar theory to either focus on the
symmetries of the critical Tokunaga family, or explore deviations from this
stiff parameterization, both of which may have physical underpinnings.

Multiple properties of the critical Tokunaga family are well justified by the
empirical evidence. We already mentioned that the coordination means that the
basin structure is determined by its Horton-Strahler order, and the Horton
prune-invariance implies that the fundamental scaling laws remain the same
after changing data resolution. Criticality ensures that the stream networks
uniformly fill the space, instead of exploding (supercritical case) or rapidly
fading off (subcritical case). The time-invariance (invariance of basin order
frequencies at different distances to the outlet) might reflect a physical
process of formation of a stream network from sources downstream, so that a
link only “knows” the information about the upstream part of the basin, yet
remains unaware of how far it is from the outlet. Deviations from this
invariance might point out to anthropogenic changes in a basin by which
various downstream alterations (dam construction, sediment aggradation, etc.)
impose upstream changes that deviate from the natural organization of a left-
alone erosional landscape. In the same vein, it would be interesting to find a
hydrogeomorphological explanation for the joint distribution of the merging
sub-basins (18).

The self-similar family extends beyond the hydrological constraints, allowing
one to study self-similar trees with edge lengths that depend on the position
within the hierarchy, arbitrary fractal dimension ${\bf d}>1$, and arbitrary
Horton exponents $R_{B}>2$ and $R_{L}>1$. For instance, the RAM might be a
suitable model for phylogenetic trees A01  or dendritic structures generated
by Diffusion Limited Aggregation (DLA). We recall that the geometric form of
the Tokunaga coefficients $T_{k}\propto c^{k}$ with $c\approx 2.72\pm 0.22$
has been known in DLA for a long time NTG97 ; O92 . It is noteworthy that the
independently estimated fractal dimension of DLA clusters, ${\bf d}=1.7\pm
0.05$ DNS86 , coincides with the fractal dimension of a critical Tokunaga tree
with $c=2.72\pm 0.22$ according to our Eq. (22): ${\bf d}=\log_{c}(2c)=1.7\pm
0.05$.

We conclude by suggesting that the understanding of the hierarchical
organization and scaling in convergent (tributary) river networks gained here
can be extended to other geomorphological processes. Important examples
include dynamic reorganization of landscapes and stream networks FGD10 ; SRF15
; WMPGC14  and scaling of peak flows MGM06 . Our results can also be extended
to study the divergent (distributary) networks of river deltas that are
commonly represented by a directed acyclic graph TLEZGRF17  – a next step in
complexity after trees examined in this work. Quantifying the structure, self-
similarity, and scaling of such graphs contributes to a still-missing unifying
theory explaining how deltaic river networks self-organize to distribute water
and sediment fluxes to the shoreline TLZF15 .

Table 1: Scaling exponents in critical Tokunaga model. Selected scaling exponents (1st column) in critical Tokunaga model expressed via the model parameter $c\geq 1$ (2nd column), fractal dimension ${\bf d}$ (3rd column), and Hack’s exponent ${\bf h}$ (4th column). Columns 5-7 show the values of the exponents in critical Tokunaga model for $c=2.3,2.4,2.5$. For comparison, column 8 shows the values estimated in the OCN model. Columns 9 summarizes estimations in the observed river networks. The agreement of the exponents of the critical Tokunaga model with $c=2.3$ (column 6) with those observed from real basins is noted.  Exponent | Expressed via | Critical Tokunaga model | OCN† | Real basins‡
---|---|---|---|---
| $c$ | ${\bf d}$ | ${\bf h}$ | $c=2.0^{*}$ | $c=2.3$ | $c=2.5$ |  |
$R_{B}=R_{M}=R_{A}$ | $2c$ | $2^{{\bf d}/({\bf d}-1)}$ | $2^{1/(1-{\bf h})}$ | 4 | 4.6 | 5.0 | 4 | 4.1 – 4.8
$R_{S}=R_{L}$ | $c$ | $2^{1/({\bf d}-1)}$ | $2^{{\bf h}/(1-{\bf h})}$ | 2 | 2.3 | 2.5 | 2 | 2.1 – 2.7
${\bf d}=\displaystyle\frac{\log R_{B}}{\log R_{L}}$ | $\log_{c}{(2c)}$ | ${\bf d}$ | ${\bf h}^{-1}$ | 2 | 1.832 | 1.756 | 2 | 1.7 – 2.0
${\bf h}=\displaystyle\frac{\log R_{L}}{\log R_{B}}$ | $\log_{2c}c$ | ${\bf d}^{-1}$ | ${\bf h}$ | 0.5 | 0.546 | 0.569 | 0.57 | 0.5 – 0.6
$\beta_{A}$ | $\log_{2c}2$ | $1-{\bf d}^{-1}$ | $1-{\bf h}$ | 0.5 | 0.454 | 0.431 | 0.43 | 0.4 – 0.5
$\beta_{\Lambda}$ | $\log_{c}2$ | ${\bf d}-1$ | ${\bf h}^{-1}-1$ | 1 | 0.832 | 0.756 | 0.8 | 0.65 – 0.9
∗ The same as the critical binary Galton-Watson model with i.i.d. exponential
edge lengths, or Shreve’s random topology model conditioned on the number of
sources.
† Average values estimated in simulated OCN basins. According to RIR01 ;
Cetal98 .
‡ According to RIR01 ; DBE94 ; Pec95 ; MRR+96 ; ZZF13 ; Rigon+96 ; TBR88 ;
TBR89 ; DR00 .

##  IV Appendices

###  IV.1 Horton-Strahler orders and Horton pruning

The importance of links and junctions in the basin hierarchy is measured by
the Horton-Strahler order $K\geq 1$ Horton45 ; Strahler . Each link and its
upstream junction have the same order. The order assignment is done in a
hierarchical fashion, from the sources downstream. Each source is assigned
order $1$. When two links of the same order $K$ merge at a junction, the
junction is assigned order $K+1$. When two links with different orders
$K_{1}>K_{2}$ merge at a junction, the largest order prevails and the junction
is assigned order $K_{1}$. The connected sequence of links and their upstream
junctions of the same order $K$ is called a stream (branch) of order $K$. We
denote by $N_{K}=N_{K}[T]$ the number of streams of order $K$ in a finite tree
$T$. The Horton-Strahler order ${\sf ord}(T)$ of a tree $T$ is the maximal
order of its links (junctions, streams). The Horton-Strahler ordering is
illustrated in Fig. 1.

The Horton-Strahler orders are closely related to the operation of Horton
pruning ${\mathcal{R}}$ that removes the source links from a basin. This
relation has been first recognized by Melton Melton59  and proved valuable in
rigorous statistical analyses of tree self-similarity Pec95 ; BWW00 ;
KZ20survey ; PG99 ; KZ16 . Formally, we consider the map
${\mathcal{R}}:{\mathcal{L}}\to{\mathcal{L}}$ that removes the source links
from a tree $T$. This may create nonbranching chains of links connected by
degree $2$ junctions – every such chain is merged into a single link. The
Horton pruning ${\mathcal{R}}$ reduces the order of each surviving stream, and
hence the basin order, by $1$. Accordingly, the order of a tree is the minimal
number of Horton prunings that completely eliminates it, as in Eq. (1). We
emphasize that the pruning cannot cut a stream in the middle – it can only
eliminate the entire stream after a finite number of iterations Melton59 .
Figure 1 illustrates the Horton pruning for the stream network of Beaver
Creek, KY – the order of this basin is $K=6$ because it is eliminated in six
Horton prunings.

###  IV.2 Asymptotic behavior of a sequence: Generating function approach

This section summarizes the basic facts about generating functions that are
the main tool in establishing asymptotic behavior of stream attributes in a
self-similar basin.

The generating function $f(z)$ of a sequence $a_{k}\geq 0$, $k=0,1,\dots$, of
non-negative real numbers is defined as a formal power series

$f(z)=\sum_{k=0}^{\infty}a_{k}z^{k},\quad z\in\mathbb{C}.$ (27)

It is known Ahlfors  that there exist such a real number $r\geq 0$ that the
series in the right hand side of (27) converges to the function $f(z)$ for any
$|z|<r$ and diverges for any $|z|>r$. The number $r$ is called the radius of
convergence of the sequence $a_{k}$; it provides notable constraints on the
asymptotic behavior of $a_{k}$. The smaller the radius of convergence, the
faster the growth of the sequence coefficients. Informally, $0<r<1$ implies
that the coefficients $a_{k}$ increase geometrically, $r>1$ that the
coefficients decrease geometrically, and $r=1$ that the coefficient vary at a
rate slower than geometric (e.g., polynomially). The values $r=0$ and
$r=\infty$ imply a faster than geometric growth or decay, respectively.

The Cauchy-Hadamar theorem Ahlfors  expresses the radius of convergence in
terms of the series coefficients:

$\frac{1}{r}=\limsup_{k\to\infty}a_{k}^{1/k}.$ (28)

Often, the radius of convergence for $a_{k}$ can be easily found from the
explicit form of $f(z)$. Specifically, if $r>0$, then the function $f(z)$ is
analytic within the disk $|z|<r$ and has at least one singularity on the
circle $|z|=r$, that is it has to diverge for at least one point on that
circle (Wilf, , Thm. 2.4.2). Thus, the radius of convergence equals to the
modulus of a singularity closest to the origin. Furthermore, recalling that
$a_{k}\geq 0$ we have

$|f(z)|=\left|\sum_{k=0}^{\infty}a_{k}\,z^{k}\right|\leq\sum_{k=1}^{\infty}a_{k}\,|z|^{k}=f(|z|),$
(29)

where the equality is only achieved for $z=|z|$. This means that the
singularity closest to the origin lies on the real axis (although there might
be other singularities with the same modulus.) This makes the search for such
a singularity much easier: one can only consider the restriction of the
function $f(z)$ on the real axis. In other words, despite the use of complex
analysis in establishing some of our results, the applied treatment of
suitable generating functions can be done in the real domain. Furthermore, if
the singularity of $f(z)$ nearest to the origin is a simple pole, then the
coefficients $a_{k}$ asymptotically form a geometric series, which we refer to
as Horton law.

######  Proposition 1 (Horton Law for a Simple Pole Sequence).

Suppose $f(z)=\sum_{i=1}^{\infty}a_{k}z^{k}$ is analytic in the disk
$|z|<\rho$ except for a single pole of multiplicity one at a positive real
value $r<\rho$. Then the sequence $a_{k}$ obeys Horton law

$\lim_{k\to\infty}a_{k}\,r^{k}=\alpha$ (30)

for some $0<\alpha<\infty$. Furthermore, if we define $g(z)=f(z)(z-r)$, then
$\alpha=-g(r)/r$.

###### Proof.

We have, for any $\Delta\in(0,r)$ Ahlfors :

$a_{k}=\frac{1}{2\pi i}\oint_{|z|=\Delta}\frac{f(z)dz}{z^{k+1}}.$ (31)

By the Residue Theorem Ahlfors , we obtain, for any $\gamma\in(r,\rho)$

$\displaystyle\frac{1}{2\pi i}\oint_{|z|=\gamma}\frac{f(z)dz}{z^{k+1}}$
$\displaystyle={\sf Res}\left(\frac{f(z)}{z^{k+1}};0\right)+{\sf
Res}\left(\frac{f(z)}{z^{k+1}};r\right)$ (32) $\displaystyle=a_{k}+{\sf
Res}\left(\frac{f(z)}{z^{k+1}};r\right).$ (33)

Therefore,

$a_{k}=\frac{1}{2\pi i}\oint_{|z|=\gamma}\frac{f(z)dz}{z^{k+1}}-{\sf
Res}\left(\frac{f(z)}{z^{k+1}};r\right),$ (34)

where

$\left|\oint_{|z|=\gamma}\frac{f(z)dz}{z^{k+1}}\right|\leq\frac{\max_{|z|=\gamma}|f(z)|}{\gamma^{k}}=o\left(r^{-k}\right).$
(35)

Consider $g(z)=(z-r)f(z)$. It is known that Ahlfors

${\sf Res}\left(f(z);r\right)=g(r),$ (36)

and hence

${\sf
Res}\left(\frac{f(z)}{z^{k+1}};r\right)=\frac{g(r)}{r^{k+1}}=\frac{g(r)}{r}r^{-k}.$
(37)

Accordingly, we obtain

$a_{k}=-\frac{g(r)}{r}r^{-k}+o(r^{-k}),$ (38)

which completes the proof. ∎

Proposition 1 is used in Sect. IV.10,IV.11 to establish Horton laws for
$\Lambda_{k}$ and $A_{k}$.

###  IV.3 Power law distribution of link attributes

Consider the value $Z_{(i)}$ of an attribute $Z$ calculated at link $i$ in a
large basin. The average number of links of order $K$ is given by
$N_{K}S_{K}$, where $S_{K}$ denotes the average number of links within a
stream of order $K$. One can heuristically approximate the frequencies of
$\\{Z_{(i)}\\}$ by using the same average value $Z_{K}$ for any link of order
$K$. Then, assuming Horton laws for $N_{K}$ and $S_{K}$ with exponents $R_{B}$
and $R_{S}$, respectively, and considering the limit of an infinitely large
basin we find

$\\#\\{i:Z_{(i)}\geq
R_{Z}^{K}\\}\approx\sum_{j=K}^{\infty}N_{j}S_{j}\propto\\!\\!\sum_{j=K}^{\infty}\\!\left(\frac{R_{S}}{R_{B}}\right)^{j}\propto\\!\left(\frac{R_{S}}{R_{B}}\right)^{K}\\!\\!\\!.$
(39)

This is a punctured (by discrete order) version of a general power law
relation of Eq. (6).

###  IV.4 Self-similar trees, Tokunaga coefficients, Horton laws

######  Definition 3 (Tokunaga coefficients).

Fix a coordinated measure $\mu$ on ${\mathcal{T}}$. For any pair $i<j$, the
Tokunaga coefficient $T_{i,j}=T_{i,j}(\mu)$ is the expected number of streams
of order $i$ per a randomly selected stream of order $j$ with respect to $\mu$
BWW00 ; Pec95 ; KZ20survey ; DR00 ; Tok78 .

We can arrange the Tokunaga coefficients for trees of a given order $K$ in an
upper triangular matrix

$\mathbb{T}_{K}=\left[\begin{array}[]{ccccc}0&T_{1,2}&T_{1,3}&\ldots&T_{1,K}\\\
0&0&T_{2,3}&\ldots&T_{2,K}\\\ 0&0&\ddots&\ddots&\vdots\\\
\vdots&\vdots&\ddots&0&T_{K-1,K}\\\ 0&0&\dots&0&0\end{array}\right].$ (40)

######  Theorem 6 (Tokunaga sequence KZ20survey ).

Suppose $\mu$ is a self-similar measure on ${\mathcal{T}}$. Then the Tokunaga
coefficients satisfy the Toeplitz property: $T_{i,i+k}=T_{k}$ for any positive
integer pair $i$, $k$. In this case the Tokunaga matrix becomes Toeplitz:

$\mathbb{T}_{K}=\left[\begin{array}[]{ccccc}0&T_{1}&T_{2}&\ldots&T_{K-1}\\\
0&0&T_{1}&\ldots&T_{K-2}\\\ 0&0&\ddots&\ddots&\vdots\\\
\vdots&\vdots&\ddots&0&T_{1}\\\ 0&0&\dots&0&0\end{array}\right].$ (41)

###### Proof.

Consider the pushforward measure ${\mathcal{R}}_{*}(\mu)$ induced on
${\mathcal{T}}$ by the Horton pruning operator:

${\mathcal{R}}_{*}(\mu)(A)=\mu\big{(}{\mathcal{R}}^{-1}(A)\big{)},\quad\forall
A\subset{\mathcal{T}}.$ (42)

Since Horton pruning decreases the order of every stream by $1$, the Tokunaga
coefficients $T_{i,j}^{R}$ computed on ${\mathcal{T}}$ with respect to
${\mathcal{R}}_{*}(\mu)$ satisfy $T_{i,j}^{R}=T_{i+1,j+1}$. The self-
similarity of $\mu$ implies $T_{i,j}^{R}=T_{i,j}$. Combining these relations,
we find $T_{i,j}=T_{i+1,j+1}$. This establishes the desired Toeplitz property
of the Tokunaga coefficients. ∎

We refer to the elements $T_{k}$ of the Tokunaga sequence as Tokunaga
coefficients, which creates no confusion with the original double-indexed
coefficients $T_{i,j}$.

######  Corollary 3 (Prune-invariance vs. Toeplitz).

Suppose $\mu$ is a coordinated measure on ${\mathcal{T}}$. Then the Toeplitz
property $T_{i,i+k}=T_{k}$ and Horton prune-invariance of Eq. (9) are
equivalent (i.e., both hold or do not hold at the same time).

###### Proof.

It has been shown in the proof of Theorem 6 that in coordinated trees both
prune-invariance and Toeplitz property take the same algebraic form
$T_{i,j}=T_{i+1,j+1}$. ∎

Consider the mean number

${\mathcal{N}}_{i}[K]={\sf E}[N_{i}(T)\,|\,{\sf ord}(T)=K]$

of streams of order $i$ in a basin of order $K$, and the mean magnitude
(number of upstream sources) $M_{i}$ of a stream of order $i$. Observe that
for a fixed $K$ the stream counts ${\mathcal{N}}_{i}[K]$ form a decreasing
sequence in $i$, and the sequence’s first term ${\mathcal{N}}_{1}[K]$
increases with $K$. At the same time, the average magnitudes $M_{i}$ form an
increasing sequence in $i$ whose first terms are independent of basin order
$K$ for any $K\geq i$. This explains the dependence on $K$ in the average
stream counts and absence of such in the average magnitudes. The definition
implies ${\mathcal{N}}_{K}[K]=M_{1}=1$ and ${\mathcal{N}}_{1}[K]=M_{K}$ for
any tree distribution. Moreover, in self-similar trees the two sequences are
deterministically related as Pec95 ; KZ16

${\mathcal{N}}_{K-j+1}[K]={\mathcal{N}}_{1}[j]=M_{k}\quad\text{for any}\quad
1\leq j\leq K.$ (43)

###  IV.5 Proof of Theorem 1

The average magnitude $M_{k}$ is the mean number of sources upstream of an
order $k$ stream. It can be represented as the sum of the magnitudes of two
order $k-1$ streams that formed this stream, plus the magnitudes of all its
side tributaries. Hence $M_{1}=1$, and

$M_{k}=2\,M_{k-1}+\sum_{i=1}^{k-1}M_{i}\,T_{k-i},\quad\text{ for }\quad k>1.$
(44)

The generating function for the average magnitudes $M_{k}$ is obtained by
multiplying both sides in (44) by $z^{k}$ and summing over $k\geq 1$:

$M(z)=\sum_{k=1}^{\infty}M_{k}\,z^{k}=z+2zM(z)+M(z)\,T(z).$

Thus,

$M(z)=\frac{z}{1-2z-T(z)}=-{z\over t(z)},$ (45)

where, according to Eq. (16) of the main text:

$t(z)=-1+2z+T(z),\quad T(z)=\sum_{i=1}^{\infty}T_{k}z^{k}.$

The function $M(z)$ is analytic with the exception of zeroes and singularities
of $t(z)$. Observe that $t(0)=-1$, and since $T_{k}\geq 0$ we have
$t(1/2)=T(1/2)\geq 0$. Furthermore, since

${d\over dz}t(z)=2+\sum_{k=1}^{\infty}\,kT_{k}\,z^{k-1}>0\,\quad\text{for all
}z\in(0,\infty),$

the equation $t(z)=0$ has a unique real root $w_{0}$ of multiplicity one in
the interval $(0,1/2]$. Let $r_{T}$ be the radius of convergence for $T(z)$,
and hence for $t(z)$. We notice that $r_{T}>w_{0}$, so the radius of
convergence of $M(z)$ coincides with the root of $t(z)$ closest to the origin.
We claim that this root is $w_{0}$. Assume otherwise, so there exists
$w\in\mathbb{C}$ such that $t(w)=0$ and $|w|<w_{0}$. Since $w_{0}$ is the
unique real root of $t(z)$ within $(0,1/2]$, $w$ must have a non-zero
imaginary part. This means that the singulatiry of $M(z)$ closest to the
origin is not on the real axis, which contradicts (29). Hence the radius of
convergence of $M(z)$ is $w_{0}$, and $w_{0}$ is a simple pole of $M(z)$.
Proposition 1 now establishes the result. ∎

###  IV.6 Exact Horton law

Assume that the Horton law for stream numbers $\mathcal{N}_{1}[K]$, and hence
for magnitudes $M_{K}$, holds exactly (recall that $M_{1}=1$):

$M_{K}=R_{M}^{K-1}.$ (46)

Then,

$M(z)=\frac{z}{1-R_{M}z}$

which leads to

$t(z)=-\frac{z}{M(z)}=-1+R_{M}z\quad\text{and}\quad T(z)=(R_{M}-2)z.$

This implies that the only self-similar model with exact Horton law
corresponds to the Tokunaga sequence

$T_{1}=R_{M}-2,\quad T_{k}=0\quad\text{for}\quad k>1.$

###  IV.7 Scalings in a critical Tokunaga tree

The critical Tokunaga model (Definition 2) has the following exact form of the
average branch counts and average magnitudes (Cor. 4 in KZ20survey ):

${\mathcal{N}}_{K-i+1}[K]=M_{i}={R_{B}^{i}+R_{B}-2\over
2(R_{B}-1)}\quad\text{for any}~{}~{}1\leq i\leq K.$ (47)

In addition, it satisfies the Horton law for the original stream counts
$N_{i}$ (KZ20survey, , Cor. 5):

${N_{i}\over N_{i+1}}\stackrel{{\scriptstyle
p}}{{\longrightarrow}}R_{B}\quad\text{ as }\,i\to\infty,$ (48)

where $\stackrel{{\scriptstyle p}}{{\longrightarrow}}$ denotes convergence in
probability BW07 . This result strengthens the statement of Theorem 1, Eq.
(II.3) that is formulated for the respective averages. Finally, the weak law
of large numbers holds for the tree order. Formally, denote by $T[K]$ a
critical Tokunaga tree of order $K$ and write $\\#T[K]$ for the number of
links in this tree. Then (KZ20survey, , Cor. 6)

${\log_{R_{B}}\\#T[K]\over K}\stackrel{{\scriptstyle
p}}{{\longrightarrow}}1\quad\text{ as }\,K\to\infty.$ (49)

Informally, this means that the tree order grows as a logarithm base $R_{B}$
of the tree size.

The identically distributed link lengths imply identically distributed local
areas, which in turn establishes the Horton law for $A_{i}$. Specifically, in
a critical Tokunaga tree we have (Appendix IV.11):

$A_{i}\sim{\rm Const.}\times M_{i},$ (50)

where $x_{i}\sim y_{i}$ stands for
$\displaystyle\lim_{i\to\infty}x_{i}/y_{i}=1$. The same approach shows that
the asymptotic of Eq. (50) holds also for the average total channel length
$L^{tot}_{i}$ upstream of a stream of order $i$, with different
proportionality constant. The asymptotic of Eq. (50) formalizes one of the key
empirical observations RIR01  that connects a physical (area $A_{i}$) and a
combinatorial (magnitude $M_{i}$) attributes of a river basin. This asymptotic
may not hold in a general self-similar (not critical Tokunaga) tree.

###  IV.8 Proofs of Corollary 1, Theorem 5, Corollary 2

###### Proof of Corollary 1.

By definition of RAM, the length of a stream of order $j$ is an exponential
random variable with rate $\lambda_{j}$. In a self-similar tree, the rate is
given by Eq. (24): $\lambda_{j}=\gamma\,\zeta^{1-j}$. This implies
$L_{j}=\lambda_{j}^{-1}=\gamma^{-1}\,\zeta^{j-1}$, which is equivalent to the
statement of Theorem 5 (Eq. (5)). ∎

###### Proof of Theorem 5.

Recall the Horton law for the average magnitude (Theorem 1, Eq. (14)) that
holds in any tree with a tamed Tokunaga sequence
($\limsup\limits_{k\to\infty}T_{k}^{1/k}<\infty$) and the Horton law for the
average length of the longest stream (Theorem 4, Eq. (23)) that holds in any
self-similar RAM tree. These laws apply to a critical Tokunaga tree of the
current statement. Furthermore, the asymptotic equivalence between the average
basin contributing area and average basin magnitude (Eq. (50)) implies the
Horton law for the average basin areas with Horton exponent $R_{M}$. Finally,
we use the general result of Eq. (5) to establish the Hack’s law (Eq. (5)) in
a critical Tokunaga tree. ∎

###### Proof of Corollary 2.

Using the definition of $\hat{t}(z)$ (Eq. (11)) and the geometric form of the
Tokunaga coefficients (Eq. (24)) we obtain $\hat{t}(z)=(1-2cz)(z-1)/(1-cz)$.
The only real root of $\hat{t}(z)$ within $(0,1/2]$ is $w_{0}=(2c)^{-1}$. By
Theorem 1, Eq. (15) we have $R_{B}=R_{M}=2c$, and Eq. (19) implies
$S_{K}=c^{K-1}$, which corresponds to $R_{S}=c$. The equality $R_{S}=R_{L}$
follows from independence of the distribution of link lengths of their
position within a basin. Finally, $R_{\Lambda}=c$ is established in Theorem 4,
Eq. (23). ∎

###  IV.9 Fractal dimension of a self-similar RAM tree

Consider a self-similar RAM tree $T$ (Theorem 3) with a Tokunaga sequence
$\\{T_{k}\\}$ satisfying $\limsup_{k\to\infty}T_{k}^{1/k}<\infty$, and
parameters $\gamma>0$ and $\zeta>1$. Below we construct a Markov tree process
$\\{\Upsilon_{K}\\}_{K=1,2,\ldots}$ corresponding to $T$ following KZ20survey
and use it to find the fractal dimension of the resulting tree in the limit of
infinite tree order. The construction below closely reproduces that of the RAM
(see the main text), but scales the edge lengths so that an infinitely large
tree has a proper fractal dimension.

Let $\Upsilon_{1}$ be an I-shaped tree of Horton-Strahler order one, with the
edge length distributed as an exponential random variable with parameter
$\gamma$. Conditioned on $\Upsilon_{K}$, the tree $\Upsilon_{K+1}$ is
constructed according to the following transition rules. We attach new leaf
edges to $\Upsilon_{K}$ at the points sampled by an inhomogeneous Poisson
point process with the intensity $\rho_{j,K}=\gamma\zeta^{K-j}T_{j}$ along the
edges of order $j\leq K$ in $\Upsilon_{K}$. We also attach a pair of new leaf
edges to each of the leaves in $\Upsilon_{K}$. The lengths of all the newly
attached leaf edges are i.i.d. exponential random variables with parameter
$\gamma\zeta^{K}$ that are independent of the combinatorial shape and the edge
lengths in $\Upsilon_{K}$. Finally, we let the tree $\Upsilon_{K+1}$ consist
of $\Upsilon_{K}$ and all the attached leaves and leaf edges.

By construction, a branch of order $j$ in $\Upsilon_{K}$ becomes a branch of
order $j+1$ in $\Upsilon_{K+1}$ after the attachment of new leave edges. The
length of order $j$ branch in $\Upsilon_{K}$ (and therefore, the length of
order $j+1$ branch in $\Upsilon_{K+1}$) is exponential random variable with
parameter $\gamma\zeta^{K-j}$. Therefore, in a tree $\Upsilon_{K+1}$, the
number $n_{1,j+1}(K+1)$ of side-branches of order one in a branch of order
$j+1$ has geometric distribution:

$\displaystyle{\sf P}\big{(}n_{1,j+1}(K+1)=r\big{)}$
$\displaystyle={\gamma\zeta^{K-j}\over\gamma\zeta^{K-j}+\rho_{j,K}}\left(\\!{\rho_{j,K}\over\gamma\zeta^{K-j}+\rho_{j,K}}\\!\right)^{r}$
$\displaystyle=\frac{1}{1+T_{j}}\left(\frac{T_{j}}{1+T_{j}}\right)^{r}$ (51)

for $r=0,1,2,\ldots$, with the mean value

${\sf
E}\big{[}n_{1,j+1}(K+1)\big{]}={\rho_{j,K}\over\gamma\zeta^{K-j}}=T_{j}.$

After $i\geq 1$ rounds of attachments the mean number $n_{i,j+i}(M)$ of side-
branches of order $i$ in a branch of order $j+i$ in a tree $\Upsilon_{M}$
(where $M=K+i$ and $K\geq j$) is

${\sf E}\big{[}n_{i,j+i}(M)\big{]}=T_{j}.$

Each tree $\Upsilon_{K}$ is distributed as a self-similar RAM tree KZ20survey
with Tokunaga sequence $\\{T_{k}\\}$ and parameters $(\gamma,\,\zeta)$,
conditioned on its Horton-Strahler order being equal to $K$, and with its edge
lengths scaled by $\zeta^{1-K}$.

Observe that by construction $\Upsilon_{K}\subset\Upsilon_{K+1}$. Accordingly,
there exists the limit space

$\Upsilon_{\infty}=\lim\limits_{K\rightarrow\infty}\Upsilon_{K}=\bigcup\limits_{K=1}^{\infty}\Upsilon_{K}.$

The self-similarity of the RAM process suggests that the limit space does not
change its statistical properties after rescaling, which corresponds here to
the Horton pruning. Let ${\bf d}$ denote its fractal dimension. That the limit
space includes at least the root branch $\Upsilon_{1}$ implies ${\bf d}\geq
1$. Assume that ${\bf d>1}$. Then, denoting the mean ${\bf d}$-dimensional
volume of $\Upsilon_{\infty}$ by ${\bf vol}$, we have

${\bf vol}=\sum\limits_{k=1}^{\infty}t_{k}{{\bf vol}\over\zeta^{{\bf d}k}}.$
(52)

This equation is obtained by splitting a tree $\Upsilon_{\infty}$ into the
subtrees attached to its highest-order branch $\Upsilon_{1}$. There is an
average of $t_{1}=T_{1}+2$ subtrees distributed as $\Upsilon_{\infty}$ scaled
by $\zeta^{-1}$. In general, for each $k$, there will be an average of $t_{k}$
subtrees distributed as $\Upsilon_{\infty}$ scaled by $\zeta^{-k}$. Scaling
the lengths by $\zeta^{-k}$ in the ${\bf d}$-dimensional space results in
scaling the volume by $\zeta^{-{\bf d}k}$. The ${\bf vol}$ term in (52) can be
cancelled out, yielding

$\hat{t}\big{(}\zeta^{-{\bf d}}\big{)}=0,$ (53)

and hence, $\zeta^{-{\bf d}}=w_{0}=R_{B}^{-1}$. This leads to (22).

###  IV.10 Horton law for $\Lambda_{k}$

If $T$ is the tree representing a stream network, then the length of the
longest stream is the height of the tree $T$, denoted by $\textsc{height}(T)$
Pitman ; KZ20survey .

Consider a tree $T$ generated by a self-similar RAM with a Tokunaga sequence
$\\{T_{k}\\}$ satisfying $\limsup\limits_{k\to\infty}T_{k}^{1/k}<\infty$, and
parameters $\gamma>0$ and $\zeta>1$. Let

$\Lambda_{k}={\sf E}\left[\textsc{height}(T)\,\Big{|}\,{\sf ord}(T)=k\right]$
(54)

that represents the mean length of the longest river stream in a basin with
the Horton-Strahler order $k$. Notice that, since
$\textsc{height}\big{(}\mathcal{R}(T)\big{)}\leq\textsc{height}(T)$,

$\displaystyle\zeta\,\Lambda_{k-1}$ $\displaystyle={\sf
E}\left[\textsc{height}\big{(}\mathcal{R}(T)\big{)}\,\Big{|}\,{\sf
ord}(T)=k\right]$ $\displaystyle\leq{\sf
E}\left[\textsc{height}(T)\,\Big{|}\,{\sf ord}(T)=k\right]=\Lambda_{k}.$ (55)

Hence, since $\Lambda_{1}=\gamma^{-1}$, we have
$\Lambda_{k}\geq\gamma^{-1}\zeta^{k-1}$. Next, let

$Y_{1},Y_{2},\ldots,Y_{N_{1}[T]}$

denote the leaf lengths in the tree $T$. Then, since

$\textsc{height}(T)\leq\textsc{height}\big{(}\mathcal{R}(T)\big{)}+\max\limits_{j=1,\ldots,N_{1}[T]}Y_{j},$

we have,

$\displaystyle\Lambda_{k}$ $\displaystyle\leq{\sf
E}\left[\textsc{height}\big{(}\mathcal{R}(T)\big{)}\,\Big{|}\,{\sf
ord}(T)=k\right]$ $\displaystyle\qquad\quad+{\sf
E}\left[\max\limits_{j=1,\ldots,N_{1}[T]}Y_{j}\,\Big{|}\,{\sf
ord}(T)=k\right]$ $\displaystyle=\zeta\,\Lambda_{k-1}\,+\gamma^{-1}{\sf
E}\left[\sum\limits_{j=1}^{N_{1}[T]}{1\over j}\,\Big{|}\,{\sf
ord}(T)=k\right]$ $\displaystyle\leq\zeta\,\Lambda_{k-1}\,+\gamma^{-1}{\sf
E}\left[1+\log\big{(}N_{1}[T]\big{)}\,\Big{|}\,{\sf ord}(T)=k\right]$
$\displaystyle\leq\zeta\,\Lambda_{k-1}\,+\gamma^{-1}+\gamma^{-1}\log\left({\sf
E}\big{[}N_{1}[T]\,\big{|}{\sf ord}(T)=k\big{]}\right)$ (56)

by Wald’s equation, the Coupon Collector Problem, and finally, the Jensen’s
inequality. Recall (Theorem 1) the Horton law for the leaf count in a self-
similar process

$\mathcal{N}_{1}[k]=M_{k}=M\,R_{B}^{k}+o\left(R_{B}^{k}\right).$

Hence, equations (IV.10) and (IV.10) imply

$0~{}\leq~{}\Lambda_{k}-\zeta\,\Lambda_{k-1}~{}\leq~{}\gamma^{-1}k\log{R_{B}}+\beta$

for some constant $\beta$, and

$0~{}\leq~{}{\Lambda_{k}\over\Lambda_{k-1}}-\zeta~{}\leq~{}\gamma^{-1}{k\log{R_{B}}+\beta\over\Lambda_{k-1}}\leq{k\log{R_{B}}+\beta\over\zeta^{k-2}}\to
0$ (57)

as $\,k\to\infty$. Accordingly,

$\displaystyle\log{\Lambda_{k}}$
$\displaystyle=\sum\limits_{j=2}^{k}\log\left({\Lambda_{k}\over\Lambda_{k-1}}\right)+\log{\Lambda_{1}}$
$\displaystyle=(k-1)\log\zeta+\sum\limits_{j=2}^{k}\log(1+\mathcal{E}_{j})-\log\gamma,$
(58)

where $0\leq\mathcal{E}_{j}\leq{(k\log{R_{B}}+\beta)\zeta^{1-k}}$, and
therefore, $\sum\limits_{j=2}^{\infty}\log(1+\mathcal{E}_{j})$ converges to a
constant. We therefore conclude that the strong Horton law holds for
$\Lambda_{k}$ with Horton exponent $R_{\Lambda}=R_{L}=\zeta$:

$\Lambda_{k}\sim\mathrm{Const.}\times\zeta^{k}.$ (59)

###  IV.11 Horton law for $A_{k}$

Assume that the mean local contributing area of a link of order $k$ equals
$\alpha_{k}$. Then the total mean contributing area $A_{K}$ of a tree of order
$K\geq 1$ is

$A_{K}=\sum\limits_{i=1}^{K}\alpha_{i}S_{i}\mathcal{N}_{i}[K],$ (60)

where $S_{k}\mathcal{N}_{k}[K]$ is the mean number of links of order $k$ in a
tree of order $K$. A convenient recursive expression is obtained by noticing
that $A_{1}=\alpha_{1}$ and

$A_{K}=2A_{K-1}+\alpha_{K}S_{K}+\sum_{i=1}^{K-1}A_{i}\,T_{K-i}\quad\text{ for
}\quad K\geq 2.$ (61)

The generating function for $A_{k}$ is given by

$A(z)=\sum\limits_{k=1}^{\infty}A_{k}z^{k}=2zA(z)+\sum\limits_{k=1}^{\infty}\alpha_{k}S_{k}z^{k}+A(z)T(z),$

which yields

$A(z)={\sum\limits_{k=1}^{\infty}\alpha_{k}S_{k}z^{k}\over
1-2z-T(z)}=-\frac{D(z)}{t(z)}=M(z)\frac{D(z)}{z}.$ (62)

Here $D(z)$ is the generating function for the mean local contributing areas
$\alpha_{k}\,S_{k}$ of streams of order $k$. Suppose that the radius of
convergence of $D(z)$ is larger than $w_{0}$. Then, by Prop. 1,

$A_{k}\sim{\rm Const.}\times w_{0}^{-k}\sim{\rm Const.}\times M_{k}.$ (63)

Consider the critical Tokunaga model. Here $\alpha_{k}=\alpha$,
$S_{k}=c^{k-1}$ and hence

$D(z)=\frac{\alpha\,z}{1-cz}$

whose radius of convergence $c^{-1}$ coincides with that of $t(z)$. Observe
that the radius of convergence of $t(z)$ must be greater than its zero, hence
$w_{0}<c^{-1}$, and so the asymptotic of (63) holds.

## References

  * (1) R. E. Horton,  Erosional development of streams and their drainage basins; hydrophysical approach to quantitative morphology.  Bulletin of Geophysical Society of America, 56, 275–370 (1945).
  * (2) I. Rodriguez-Iturbe and A. Rinaldo,  Fractal river basins: chance and self-organization.  Cambridge University Press (2001).
  * (3) J. W. Kirchner,  Statistical inevitability of Horton’s laws and the apparent randomness of stream channel networks.  Geology, 21(7) 591–594 (1993).
  * (4) R. L. Shreve ,  Statistical law of stream numbers.  J. Geol., 74(1) 17–37 (1966).
  * (5) G. Burd, E. C. Waymire, and R. D. Winn,  A self-similar invariance of critical binary Galton-Watson trees.  Bernoulli Society for Mathematical Statistics and Probability, 6(1), 1–21 (2000).
  * (6) J. Pitman,  Combinatorial Stochastic Processes.  Ecole d’été de probabilités de Saint-Flour XXXII-2002.  Lectures on Probability Theory and Statistics, Springer, (2006).
  * (7) H. De Vries, T. Becker, and B. Eckhardt,  Power law distribution of discharge in ideal networks.  Water Resources Research, 30(12), 3541–3543 (1994).
  * (8) S. D. Peckham,  New Results for Self-Similar Trees with Applications to River Networks.  Water Resources Research, 31(1), 1023–1029 (1995).
  * (9) A. E. Scheidegger,  A stochastic model for drainage patterns into an intramontane treinch.  Hydrological Sciences Journal, 12(1), 15–20 (1967).
  * (10) H. Takayasu, I. Nishikawa, and H. Tasaki,  Power-law mass distribution of aggregation systems with injection.  Physical Review A, 37(8), 3110 (1988).
  * (11) A. Rinaldo, I. Rodriguez-Iturbe, R. Rigon, R. L. Bras, E. Ijjasz-Vasquez, and A. Marani,  Minimum energy and fractal structures of drainage networks.  Water Resources Research, 28(9), 2183–2195 (1992).
  * (12) R. Rigon, A. Rinaldo, I. Rodriguez-Iturbe, R. L. Bras, and E. Ijjasz-Vasquez,  Optimal channel networks: a framework for the study of river basin morphology.  Water Resources Research, 29(6), 1635–1646 (1993).
  * (13) A. Rinaldo, I. Rodriguez-Iturbe, R. Rigon, E. Ijjasz-Vasquez, and R. L. Bras,  Self-organized fractal river networks.  Physical Review Letters, 70(6), 822 (1993).
  * (14) A. Maritan, A. Rinaldo, R. Rigon, A. Giacometti,and I. Rodriguez-Iturbe,  Scaling laws for river networks.  Physical Review E, 53(2), 1510 (1996).
  * (15) A. Rinaldo, R. Rigon, J. R. Banavar, A. Maritan, I. Rodriguez-Iturbe,  Evolution and selection of river networks: Statics, dynamics, and complexity.  Proceedings of the National Academy of Sciences, 111(7), 2417–2424 (2014).
  * (16) P. Balister, J. Balogh, E. Bertuzzo, B. Bollobás, G. Caldarelli, A. Maritan, R. Mastrandrea, R. Morris, and A. Rinaldo,  River landscapes and optimal channel networks.  Proceedings of the National Academy of Sciences, 115(26), 6548–6553 (2018).
  * (17) V. K. Gupta and E. C. Waymire,  Statistical self-similarity in river networks parameterized by elevation.  Water Resources Research, 25(3) 463–476 (1989).
  * (18) V. K. Gupta, B. M. Troutman, and D. R. Dawdy,  Towards a nonlinear geophysical theory of floods in river networks: an overview of 20 years of progress.  Nonlinear Dynamics in Geosciences (pp. 121–151), Springer, New York, NY (2007).
  * (19) D. L. Turcotte,  Fractals and chaos in geology and geophysics.  Cambridge University Press, (1997).
  * (20) W. I. Newman, D. L. Turcotte and A. M. Gabrielov,  Fractal trees with side-branching.  Fractals, 5, 603–614 (1997).
  * (21) S. D. Peckham and V. K. Gupta,  A reformulation of Horton’s laws for large river networks in terms of statistical self-similarity.  Water Resources Research, 35(9), 2763–2777 (1999).
  * (22) P. Flajolet, J.-C. Raoult, and J. Vuillemin,  The number of registers required for evaluating arithmetic expressions.  Theoretical Computer Science 9(1) 99–125 (1979).
  * (23) M. Drmota and H. Prodinger,  The register function for t-ary trees ACM.  Transactions on Algorithms 2 (3) 318–334 (2006).
  * (24) M. Baiesi and M. Paczuski,  Scale-free networks of earthquakes and aftershocks.  Physical Review E, 69(6), 066106 (2004).
  * (25) J. R. Holliday, D. L. Turcotte, and J. B. Rundle,  Self-similar branching of aftershock sequences.  Physica A: Statistical Mechanics and its Applications, 387(4) 933–943 (2008).
  * (26) M. R. Yoder, J. Van Aalsburg, D. L. Turcotte, S. G. Abaimov, and J. B. Rundle,  Statistical variability and Tokunaga branching of aftershock sequences utilizing BASS model simulations.  Pure and Applied Geophysics, (2013) 170(1-2) 155–171.
  * (27) I. Zaliapin, and Y. Ben-Zion,  Earthquake clusters in southern California I: Identification and stability.  Journal of Geophysical Research: Solid Earth, 118(6), 2847–2864 (2013).
  * (28) G. S. Kassab,  The coronary vasculature and its reconstruction.  Annals of Biomedical Engineering, 28(8) 903–915 (2000).
  * (29) F. Cassot, F. Lauwers, C. Fouard, S. Prohaska, and V. Lauwers-Cances,  A novel three-dimensional computer-assisted method for a quantitative study of microvascular networks of the human cerebral cortex.  Microcirculation, 13(1) 1–18 (2006).
  * (30) E. H. Campbell Grant, W. H. Lowe, and W. F. Fagan,  Living in the branches: population dynamics and ecological processes in dendritic networks.  Ecology Letters, 10(2) 165–175 (2007).
  * (31) D. L. Turcotte, J. D. Pelletier, and W. I. Newman,  Networks with side-branching in biology.  Journal of Theoretical Biology, 193(4), 577–592 (1998).
  * (32) Y. Kovchegov and I. Zaliapin,  Random Self-Similar Trees: A mathematical theory of Horton laws.  Probability Surveys, 17, 1–213 (2020).
  * (33) A. N. Strahler,  Quantitative analysis of watershed geomorphology.  Trans. Am. Geophys. Un., 38 913–920 (1957).
  * (34) M. A. Melton,  A derivation of Strahler’s channel-ordering system.  The Journal of Geology, 67(3), 345–346 (1959).
  * (35) L. Leopold, M. G. Wolman, and J. Miller,  Fluvial processes in geomorphology.  Dover Publications, Inc., New York, USA (1992).
  * (36) D. G. Tarboton,  Fractal river networks, Horton’s laws and Tokunaga cyclicity.  Journal of hydrology, it 187(1) 105–117 (1996).
  * (37) V. K. Gupta and E. C. Waymire,  Some mathematical aspects of rainfall, landforms and floods.  In O. E. Barndorff-Nielsen, V. K. Gupta, V. Perez-Abreu, E. C. Waymire (eds) Rainfall, Landforms and Floods, Singapore: World Scientific, Singapore (1998).
  * (38) S. Zanardo, I. Zaliapin, and E. Foufoula-Georgiou,  Are American rivers Tokunaga self-similar? (in New results on fluvial network topology and its climatic dependence).  Journal of Geophysical Research: Earth Surface, 118, 1–18 (2013).
  * (39) O. J. Mesa,  Cuatro modelos de redes de drenaje.  Revista de la Academia Colombiana de Ciencias Exactas, Físicas y Naturales, 42(165), 379–391 (2018).
  * (40) J. T. Hack,  Studies of longitudinal stream profiles in Virginia and Maryland.  US Government Printing Office, 294 (1957).
  * (41) R. Rigon, I. Rodriguez-Iturbe, A. Maritan, A. Giacometti, D. Tarboton, and A. Rinaldo,  On Hack’s Law.  Water Resources Research, 32(11), 3367–3374 (1996).
  * (42) D. G. Tarboton, R. L. Bras, and I. Rodriguez-Iturbe,  The fractal nature of river networks.  Water Resources Research, 24, 1317–1322 (1988).
  * (43) I. Rodriguez-Iturbe, E. J. Ijjasz-Vasquez, R. L. Bras, and D. G. Tarboton,  Power law distributions of discharge mass and energy in river basins.  Water Resources Research, 28(4) 1089–1093 (1992).
  * (44) Y. Kovchegov and I. Zaliapin,  Horton Law in Self-Similar Trees.  Fractals, 24, 1650017 (2016).
  * (45) P. La Barbera and R. Rosso,  On the fractal dimension of stream networks.  Water Resources Research, 25(4), 735–741 (1989).
  * (46) D. G. Tarboton, R. L. Bras, and I. Rodriguez-Iturbe,  Scaling and elevation in river networks.  Water Resources Research, 25(9), 2037–2051 (1989).
  * (47) E. Tokunaga,  Consideration on the composition of drainage networks and their evolution.  Geographical Reports of Tokyo Metropolitan University, 13, 1–27 (1978).
  * (48) S. A. Veitzer and V. K. Gupta,  Random self-similar river networks and derivations of generalized Horton Laws in terms of statistical simple scaling.  Water Resour. Res., 36(4) 1033–1048 (2000).
  * (49) Y. Kovchegov and I. Zaliapin,  Tokunaga self-similarity arises naturally from time invariance.  Chaos: An Interdisciplinary Journal of Nonlinear Science, 28(4) 041102 (2018).
  * (50) E. V. Chunikhina,  Entropy rates for Horton self-similar trees.  Chaos, 28(8), 081104 (2018).
  * (51) E. V. Chunikhina,  Information theoretical analysis of self-similar trees.  Ph.D. thesis, Oregon State University (2018).
  * (52) P. Dodds and D. Rothman,  Scaling, Universality, and Geomorphology.  Ann. Rev. Earth and Planet. Sci., 28 571–610 (2000).
  * (53) V. K. Gupta, S. L. Castro, and T. M. Over,  On scaling exponents of spatial peak flows from rainfall and river network geometry.  Journal of Hydrology, 187(1-2) 81–104 (1996).
  * (54) M. Cieplak, A. Giacometti, A. Maritan, A. Rinaldo, I. Rodriguez-Iturbe, and J. R. Banavar,  Models of fractal river basins.  Journal of Statistical Physics, 91(1–2) 1–15 (1998).
  * (55) J. D. Pelletier and D. L. Turcotte,  Shapes of river networks and leaves: are they statistically similar?  Philosophical Transactions of the Royal Society of London B: Biological Sciences, 355(1394), 307–311 (2000).
  * (56) B. Lashermes and E. Foufoula-Georgiou,  Area and width functions of river networks: New results on multifractal properties.  Water Resour. Res., 43, W09405, (2007).
  * (57) V. K. Gupta and E. C. Waymire,  Multiscaling properties of spatial rainfall and river flow distributions.  Journal of Geophysical Research: Atmospheres, 95(D3), 1999–2009 (1990).
  * (58) H. W. Watson, F. Galton,  On the probability of the extinction of families.  The Journal of the Anthropological Institute of Great Britain and Ireland, 4, 138–144 (1875).
  * (59) D. J. Aldous,  Stochastic models and descriptive statistics for phylogenetic trees, from Yule to today.  Statistical Science, 16(1), 23–34 (2001).
  * (60) P. Ossadnik,  Branch order and ramification analysis of large diffusion-limited-aggregation clusters,  Physical Review A, 45(2), 1058 (1992).
  * (61) G. Daccord, J. Nittmann, and H. E. Stanley,  Radial viscous fingers and diffusion-limited aggregation: Fractal dimension and growth sites.  Physical Review Letters, 56(4), 336 (1986).
  * (62) E. Foufoula-Georgiou, V. Ganti, and W. E. Dietrich,  A non-local theory for sediment transport on hillslopes.  J. Geophys. Res. - Earth Surface, 115, F00A16 (2010).
  * (63) A. Singh, L. Reinhardt, and E. Foufoula-Georgiou,  Landscape reorganization under changing climatic forcing: results from an experimental landscape.  Water Resour. Res., 51(6), 4320–4337 (2015).
  * (64) S. D. Willett, S. W. McCoy, J. T. Perron, L. Goren, and C. Y. Chen,  Dynamic reorganization of river basins.  Science, 343(6175) (2014).
  * (65) R. Mantilla, V. K. Gupta, and O. J. Mesa,  Role of coupled flow dynamics and real network structures on Hortonian scaling of peak flows.  Journal of Hydrology, 322(1-4), 155–167 (2006).
  * (66) A. Tejedor, A. Longjas, D. A. Edmonds, I. Zaliapin, T. T. Georgiou, A. Rinaldo, and E. Foufoula-Georgiou,  Entropy and optimality in river deltas.  Proceedings of the National Academy of Sciences, 114(44), 11651–11656 (2017).
  * (67) A. Tejedor, A. Longjas, I. Zaliapin, and E. Foufoula-Georgiou,  Delta channel networks: 1. A graph-theoretic approach for studying connectivity and steady state transport on deltaic surfaces.  Water Resources Research, 51(6), 3998–4018 (2015).
  * (68) L. V. Ahlfors,  Complex analysis: an introduction to the theory of analytic functions of one complex variable.  New York, London, 177 (1953).
  * (69) R. N. Bhattacharya and E. C. Waymire,  A basic course in probability theory (Vol. 69), Springer, New York (2007).
  * (70) H. S. Wilf,  Generatingfunctionology.  Philadelphia, PA, USA (1992).