---
abstract: 'The dynamics of a 2[*D*]{} site percolation model on a square lattice is studied using the hierarchical approach introduced by Gabrielov et al., [*Phys. Rev. E*]{}, [**60**]{}, 5293-5300, 1999. The key elements of the approach are the tree representation of clusters and their coalescence, and the Horton-Strahler scheme for cluster ranking. Accordingly, the evolution of percolation model is considered as a hierarchical inverse cascade of cluster aggregation. A three-exponent time-dependent scaling for the cluster rank distribution is derived using the Tokunaga branching constraint and classical results on percolation in terms of cluster masses. Deviations from the pure scaling are described. An empirical constraint on the dynamics of a rank population is reported based on numerical simulations.'
author:
- 'Ilya Zaliapin, Henry Wong, Andrei Gabrielov'
title: 'Inverse cascade in percolation model: hierarchical description of time-dependent scaling'
---

Introduction
============

Percolation model is probably the simplest and best studied system that experiences (geometrical) phase transition of the second kind [@SA]. It is widely used as a toy model for spatially distributed stochastic processes, such as diffusion in disordered media, forest fires, gelation, semiconduction, etc. [@Sornette04; @SA]. Importantly for our study, percolation model presents a transparent mechanism of the process of hierarchical aggregation (coagulation). This process has been actively employed for describing the essential properties of material fracture and earthquake nucleation [@BSLA97; @BSL97; @GKZN00; @NS90; @NG91; @NTG95; @SB01; @ZKG03], starting from the pioneering works of Allegre et al. [@ALP82] and Newman and Knopoff [@NK82; @NK83; @KN83]. In this paper we describe the evolution of percolation model in terms of an inverse cascade of hierarchical cluster aggregation.

An early idea of hierarchical aggregation was introduced by Newman and Knopoff in the “crack-fusion” model for repetitive cycles of large earthquakes [@NK82; @NK83; @KN83; @NK90]. Their model focused on processes of small cracks fusions into successively larger ones, accommodating the influence of mainshocks and aftershocks, juvenile crack genesis from tectonic stresses, crack healing, and anelastic-creep induced time delays, plus other effects. Turcotte et al. [@Tur+99] have reinstated this line of research considering a log-binned description of hierarchical aggregation and performing numerical tests to study its scaling properties. Gabrielov et al. [@GNT99] first have employed the Horton-Strahler hierarchical ranking [@Horton45; @NTG97] to construct an exactly solvable model of a general inverse cascade process. The Horton-Strahler ranks (see Sect. \[HS\]) that came from hydrology and have been not well known in physical applications happened to be more natural than cluster masses (sizes, areas) in describing the aggregation process. Moreover, the ranks are shown essential for formulating the analytical models [@GNT99]. Recent efforts deal with studying the aggregation dynamics and its various scalings via exactly solvable hierarchical models and extensive simulations [@MNTG04].

Below we focus on the evolution of the first spanning cluster in the the classical site-percolation model, and decribe it as a consecutive hierarchical fusion of smaller clusters into larger ones. Noteworthy, we are interested not in a final solution of a percolation state, but in an evolutionary path leading from the juvenile single-particle clusters to a self-similar population of clusters of arbitrary large size (limited by the finiteness of the lattice), the percolation cluster included. Thus we depart from the steady-state assumption of [@GNT99; @MNTG04; @Tur+99] as well as from the asymptotic focus on the percolation onset typical for the classical percolation studies [@SA].

Specifically, we follow [@GNT99] and represent each cluster by a time-oriented tree that reflects the history of cluster formation. The model dynamics is then described in terms of the corresponding trees using the well-developed theory of hierarchical scaling complexities [@BP97; @NTG97]. An important role is played by the the Horton-Strahler scheme that provides a natural ranking for the tree-based structures. Another important element is the Tokunaga classification that defines a special subclass of trees with self-similar branching. A large number of hierarchies observed in nature are shown to be Tokunaga trees [@NTG97]; this is also the case for the clusters in percolation model [@GNT99; @MNTG04]. We use the Tokunaga constraint together with classical results on percolation dynamics (in terms of cluster masses) to derive time-dependent scaling laws for rank distribution of clusters. Importantly, we report a three-exponent scaling for the dynamics of a population of clusters of a given rank, in deviation from the two-exponent scaling well-known for the population of a given mass [@SA; @Margolina+84]. We also analyze deviations from the pure scaling and confirm our results by numerical simulations.

The inverse cascades and aggregation (coagulation) processes are important for evolution of many natural hazardous processes: earthquakes, landslides, and forest fires are argued to follow the hierarchical aggregation dynamics [@TMGR02; @MNTG04]. A general review of the theory and models of kinetics of irreversible aggregation is given by Leyvraz [@Leyvraz03]. An alternative approach to analytical modeling, based on ideas from [@GNT99], but using equations that are consistent with the mass action law of chemical kinetics, can be found in da Costa et al. [@CGW02].

The paper is organized as follows. The percolation model is described in Sect. \[model\]; this section also introduces tree representation of clusters and the Horton-Strahler ranking. In Sect. \[MRD\] we derive the average mass of clusters of a given rank using the Tokunaga constraint on cluster branching. This result will be actively used in consecutive sections. Section \[RD\] is devoted to the time-dependent rank distribution of clusters. First (Sect. \[rankdist\]), we establish the exponential rank distribution at percolation using the result of Sect. \[MRD\]. We then proceed with time-dependent rank distribution; Sect. \[dynrankdist\] introduces the three-exponent scaling for ranks and compares it to the well-known Stauffer’s two-exponent scaling for cluster masses. Scaling for ranks averaged over the entire evolution of the percolation cluster is derived in Sect. \[AS\]; this result is motivated by the heuristic studies that typically use averaged observations on a system. Time-dependent finite-size corrections to the established scalings are described in Sect. \[corrections\]. Our study of rank distributions is concluded in Sect. \[MDGR\] by describig the time-dependent behavior of the total mass of clusters of a given rank. Sect. \[CFS\] analyzes fractal properties of clusters and reports sharp increase of cluster fractal dimension in the vicinity of percolation. Sect. \[DC\] uses simulations to establish a notable constraint on the dynamics of a rank populations. The results are discussed in Sect. \[discussion\].

Model
=====

Dynamics
--------

We consider the classical 2[*D*]{} site-percolation model [@SA]. The model dynamics starts with an empty $L\times L$ square lattice. At each step a particle is dropped into a randomly chosen unoccupied site; thus each site can be either occupied by one and only one particle or empty. Two sites are considered [*neighbors*]{} if they share one side; each site on a square lattice has four neighbors. Cluster is defined as a group of occupied neighbor sites [@SA]. Time refers to the steps at which particles drop onto the lattice. Since we do not have annihilation of particles, time is formally equivalent to the number of particles on the lattice. It is convenient to normalize time by the lattice size $L^2$ so it varies from $\rho=0$ at the start to $\rho=1$ when all sites are occupied. During the system evolution, occupied sites start to aggregate and clusters begin to form. Once a sufficient number of particles is accumulated, a percolation cluster is formed connecting the opposite sides of the lattice vertically and/or horizontally.

The density $\rho$ increases monotonically from zero to its critical value $\rho_c$ at percolation. For an infinite lattice $\rho_c \approx 0.59274606$ [@NZ01], while for a finite lattice it is smaller [@SA]: $$\label{rhofin}
\rho_c(L) = \rho_c - cL^{-3}.$$

Many phenomena encountered in the percolation model mimic what we see when the phase transitions of the second kind occur. Note however that these phenomena are of purely geometrical and statistical rather than physical nature. Indeed, the physical percolation theory is largerly predicated in this geometrical model and there are many empirical links between them; this is why the percolation model is said to be an example of the [*geometrical*]{} phase transition of the second kind, and why its nomenclature emerges from that of the physical critical phenomena.

The theoretical description of the percolation dynamics is conventionally given in terms of the cluster masses [@SA]; and most of the universal scalings – a benchmark of phase transitions of second kind – deal with parameters expressed via the mass distribution of clusters. However, if one is interested in analytical description of the aggregation process, the mass description happens to be inferior to the hierarchical rank approach [@GNT99; @MNTG04]. Properly defined ranks not only allow one to construct exactly solvable models of aggregation, but also they are more feasible for observations in practice. In addition, they reflect the individual history of cluster formation. Below we follow the hierarchical approach of Gabrielov et al. [@GNT99] to study the percolation dynamics.

Tree representation of clusters
-------------------------------

Each cluster in our model is represented by a tree that reflects the time-dependent formation of a cluster (its history), and is a subject for quantitative analysis. Specifically, each one-particle cluster is represented by a trivial tree consisting of a single node. When two clusters are merged together their trees are also merged by adding a new node (parent) for which they become children (and siblings to each other.) In our model, the coalescence of two or more clusters can only be materialized by adding to the lattice a new particle which will be a neighbor to one or more existing clusters. Figure \[fig\_model1\]a illustrates the four possible types of coalescence. We call $k$-coalescence in a situation when a newly dropped particle (marked [**N**]{} in the figure) is a neighbor to $k$ existing clusters (gray numbered sites). Numerical simulations on a square lattice with $L=2,000$ suggest the following relative frequencies $Q_k$ of $k$-coalescences: $Q_1\approx0.628$, $Q_2\approx0.318$, $Q_3\approx0.052$, $Q_4\approx0.002$. Figures \[fig\_model1\]b,c illustrate how a tree is formed for different coalescence types. There are two basic situations: When a new particle is a neighbor to only one existing cluster, it is considered as an individual one-particle cluster that is connected to the existing one. The connecting node of the tree in this case does not correspond to a particle on the lattice (panel b). When a new particle is droped in a neighbor to two, three, or four existing clusters, it is not condidered as an individual cluster. Instead, it corresponds to the connecting node in the tree (panel c). Thus, the connecting node in a tree may or may not correspond to a lattice particle depending on the coalescence type. The branching parameter (number of children for a given parent) of a tree for any cluster varies between 2 and 4. Note that both 1- and 2-coalescences result in merging only two clusters; accordingly, most of the observed coalescences (about 95%) involve only two clusters while coalescence of three or four clusters is extremely rare.

The consecutive process of tree formation for a simple four-particle cluster is illustrated in Fig. \[fig\_model2\]. Importantly, the individual evolution of a cluster is crucial in constructing the corresponding hierarchical tree. To construct the tree one needs to consider all consecutive coalescences that have formed the cluster, not only its final shape. Therefore, it is clear that the same tree may correspond to clusters of different shape: Figure \[fig\_model3\]a shows two 11-particle clusters that both correspond to the same tree shown in panel b. Therefore, working with trees, we unavoidably narrow the information about the cluster population. Notice however both trees capture an excessively larger amount of information than mere cluster masses. Summing up, the time evolution of a cluster is neccesary and sufficient to uniquely determine the corresponding tree, while the inverse is not true. The problems of describing the set of trees that might correspond to a given cluster, and the set of clusters that correspond to a given tree is beyond the scope of this paper. Next, we describe the ranking of clusters, presenting a conventional alternative to the logarithmic binning of cluster masses.

Horton-Strahler ranking {#HS}
-----------------------

The appropriate ordering of trees (clusters) is very important for meaningful description and analysis of the model dynamics. The problem of such an ordering is not trivial since the clusters may grow and coalesce in a variety of peculiar ways. An advantageous way to solve this problem is given by the Horton-Strahler topological classification of ramified patterns [@Horton45; @Strahler; @BP97] illustrated in Fig. \[fig\_model3\]b: One assigns ranks to the nodes of a tree, starting from $r=1$ at leaves (clusters consisting of one particle.) When two or more clusters with ranks $r_i$, $i=1,\dots,n$ merge together, a new cluster is formed with the rank [@BP97]: $$r=\left\{\begin{array}{ll}
r_1+1,&\mbox{ if }r_i=r_1~\forall~i=1,\dots,n\\
\max{(r_i)},&\mbox{ otherwise.}
\end{array}\right.$$ The rank of a cluster is that of the root of the corresponding tree. It is possible to consider an alternative definition of ranks: When at least two clusters with rank $r$ coalesce, and other coalescing clusters have a lower rank, the rank of a new cluster becomes $r+1$. Clearly, the two definition coincide when only two clusters coalesce. The results reported in this paper are independent of the particular definition, since coalescence of more than two clusters (especially of high ranks) is a rare event.

Originally introduced in geomorphology by Horton [@Horton45] and later refined by Strahler [@Strahler], this classification is shown to be inherent in various geophysical, biological, and computational applications [@BP97; @GNT99; @MNTG04; @NTG97; @Toro01; @TPN98].

Mass-rank distribution {#MRD}
======================

In this section we derive the distribution of the average mass $m_r$ of rank $r$ clusters. It will be used consequently to connect various mass and rank scaling laws. First, we define the branching ratio $T_{ij}$ for a given cluster (tree) as the number $N_{ij}$ of subclusters (nodes) of rank $i$ that joined subcluster (node) of rank $j$, averaged over subclusters (nodes) of rank $j$ [@NTG97; @Tokunaga]: $$T_{ij}=\frac{N_{ij}}{N_j}.$$

Next we note that the mass of a rank $r$ cluster is the sum of two $r-1$ cluster masses that formed the cluster (we ignore the possibility for three or more clusters to coalesce at the same step), plus a unit mass of a joining particle, plus the mass of all the lower-rank clusters that joined the considered cluster, hence: $$\begin{aligned}
\label{mr}
m_1&=&1\nonumber\\
m_2&=&(2\,m_1+P)+T_{12}(m_1+P)\nonumber\\
m_3&=&(2\,m_2+1)+T_{23}(m_2+1)+T_{13}(m_1+P)\nonumber\\
&\dots&\nonumber\\
m_k&=&(2\,m_{k-1}+1)+\sum_{i=1}^{k-1}T_{k-i\,k}(m_{k-i}+1)
-(1-P)T_{1\,k},~k\ge 3.\end{aligned}$$ Here the coefficient $P$ addresses the possibility for a one-particle cluster to join another cluster in two ways: via a one-particle connector (with probability $P$) or directly (with probability $1-P$); the clusters with $r>2$ can only join other clusters using a one-particle connector.

It was predicted by Gabrielov et al. [@GNT99] and later confirmed by simulations [@MNTG04] that clusters in percolation model obey the Tokunaga scaling [@Tokunaga] asymptotically in $k$: $$\label{tokunaga}
T_{i\,i+k}=T_k=s_0s^{k-1}.$$ This rewrites Eq. (\[mr\]) for $k\ge 3$ as $$m_k=(2\,m_{k-1}+1)+\sum_{i=1}^{k-1}T_{i}(m_{k-i}+1)
-(1-P)T_{k-1}.$$ Assuming the mass-rank relation in the form $m_r=c^{r-1}$, $c>1$ we obtain $$\begin{aligned}
\label{mr1}
c\,^{k-1}&=&2c\,^{k-1}+1+\sum_{i=1}^{k-1}s_0s\,^{i-1}
\left(c\,^{k-i-1}+1\right)-
(1-P)s_0\,s\,^{k-1}\nonumber\\
&=&c\,^{k-2}\left[
2+\frac{1}{c\,^{k-2}}
+s_0\frac{1-\left(s/c\right)^{k-1}}
{1-s/c}
+\frac{s_0}{c\,^{k-2}}\,\frac{s\,^{k-1}-1}{s-1}
-(1-P)s_0\left(s/c\right)^{k-2}\right].\nonumber\end{aligned}$$ It is easily checked that this equation has a solution only if $c>s$; thus $s/c<1$ and for large $k$ then follows $$c^{k-1} = c^{k-2}\left[2+\frac{s_0}{1-s/c}\right]$$ leading to the final equation $$c^2-c(2+s+s_0)+2s=0$$ with solution: $$\label{c}
c=\frac{2+s+s_0\pm\sqrt{(2+s+s_0)^2-8s}}{2}.$$

Remarkably, the model of Gabrielov et al. [@GNT99] predicts in a Euclidean (assuming clusters of regular, non-fractal, shape) limit of an inverse cascade model $$s_0\approx 0.55495813,~~ s=1/s_0\approx1.80193774,
~~{\rm and~~} c=1/s_0^2\approx3.24697602.$$ The Eq. (\[c\]) in this case gives $c(s_0,s)=3.24697960$ (this is the only solution such that $c>s$), which is remarkably close (6 digits) to the result of [@GNT99]. Furthermore, the non-Eucledian (assuming fractal shape of clusters) steady-state simulations of Morein et al. [@MNTG04] suggest $$s\approx3.0253, \qquad s_0\approx0.6993,
\qquad c\approx4.325,$$ which exactly solves Eq. (\[c\]). We found it quite amazing that our complimentary set of assumptions used to derive (\[c\]) lead to the same numerical results as analytical study [@GNT99] and simulations of [@MNTG04]. This suggests an underlying connection between our approaches to describe the hierarchical aggregation.

The observed mass-rank distribution of clusters at percolation is shown in Fig. \[fig\_MR\]; it obeys the exponential relation $$\label{MR}
m_r= 10^{\,\gamma (r-1)} = c^{\,r-1},$$ with $\gamma\approx 0.625$, $c=10^{\gamma}\approx 4.2$. Our simulation suggest that the mass distribution within a given rank is approximately lognormal (not shown) with the mean given by Eq. (\[MR\]) and a rank-independent standard deviation.

The relation (\[MR\]) is a key element in our further analysis. As we will show, the distribution of cluster ranks at percolation (Sect. \[rankdist\]) and its finite-size corrections (Sect. \[corrections\]) are obtained from the corresponding classical laws for masses by simple substituting the relation (\[MR\]). At the same time, one of the most important results: the time dependent rank distribution can not be obtained this way and requires an additional treatment (Sect. \[dynrankdist\]).

The exponential relation of Eq. (\[MR\]) happens to be valid over the entire time interval $0<\rho\le\rho_c$. The corresponding dynamics of $c(\rho)$ is shown in Fig. \[fig\_c\]: it grows with time from about 2.0 at the earliest stages to 4.2 at percolation. This growth reflects the fact that clusters become more weighty with time due to coupling with the clusters of lower ranks (which does not change the rank but increases the mass). The growth is not monotonous; it is accompanied by pronounced log-periodic oscillations which are associated with creation of new ranks. The log-periodic oscillations that accompany general power-law increase of observed parameters have been found in many systems including hierarchical models of defect development [@NTG95], biased diffusion on random lattices [@SS98], diffusion-limited aggregation (DLA) [@Sor+96], and others. Log-periodic oscillations can be naturally explained by the Discrete Scale Invariance (DSI) [@Sornette04], which occurs in a system whose observables scale only for a discrete set of values. A famous example of DSI is given by the Cantor set that pocesses a discrete scale symmetry: In order to superimpose its scaled image onto the original, one has to stretch it by the discrete factors $3^n$, $n=1,2,\dots$, not a continuous set of values. The Cantor set and percolation belong to systems with built-in geometrical hierarchy, leading to DSI. In our particular system, ranks take only a countable set of values. Creation of new ranks necessarily disrupt the system in a discontinuous way resulting in the log-periodicity.

Now we return to the numerical value of parameter $c$. In steady-state simulations of [@MNTG04] $c = 4.325$, which is reasonably close to what we observe at percolation. Recall that the models of [@GNT99; @MNTG04] use the “fractal correction” $\epsilon$ to the cluster shape; this correction affects the rate $r_{ij}$ of clusters coalescence: $$r_{ij}\approx \epsilon^{-|j-i|}L_iL_j,$$ where $L_i$ is the total boundary size of the clusters of rank $i$. The correction $\epsilon$ can be expressed as $$\epsilon=\frac{1}{\sqrt{c}}\,\frac{c-1}{c-2},$$ which, together with results of Fig. \[fig\_c\], shows that in the percolation model $\epsilon$ decreases in time passing the Euclidean limit $\epsilon=1$ [@GNT99] at $(\rho_c-\rho)\approx 0.14$ and approaching the steady-state “fractal” $\epsilon\approx 0.68$ [@MNTG04] at $\rho=\rho_c$. The interval $2< c \le 4.2$ observed during $0<\rho\le \rho_c$ corresponds to $0.68\le\epsilon<\infty$.

Rank distribution {#RD}
=================

This section is devoted to establishing various time-dependent scaling laws for clusters of a given rank. We will see that it is tipically impossible to derive such laws by applying the mass-rank relation (\[MR\]) to the coresponding well-known laws for cluster masses. This illustrates an original character and richness of the rank description and prompts for developing new methods of analysis. We start with the simplest problem: rank distribution at percolation.

Distribution at percolation {#rankdist}
---------------------------

We start recalling the well-known cluster mass distribution at percolation [@SA]: $$\label{GRm}
n_m(\rho_c)\sim q_0\,m^{-\tau},$$ where $n_m(\rho_c)$ is the number of clusters of mass $m$ per lattice site, and the Fisher exponent $\tau=187/91\approx 2.05$ is universal for 2[*D*]{} systems [@Fisher67; @SA]. Figure \[fig\_GRm\] illustrates the mass distribution at percolation for a system with $L=2000$; to smooth out statistical fluctuations it shows the number of clusters with mass equal to or larger than $m$: $\sum_{m'\ge m} n_m(\rho_c)$. Equation (\[GRm\]) suggests the slope $\tau-1\approx\,1.05$, while the observed slope $0.96$ is somewhat less than that. This is due to the impact of two concurrent phenomena: so-called “deviation from scaling” at small $m$ [@Hoshen+79] and finite-size effects at large $m$ [@Margolina+84; @Hoshen+79]; they are discussed below in Sect. \[corrections\].

Now, we use Eq. (\[GRm\]) to derive the distribution of the number $n_r(\rho_c)$ of the clusters of rank $r$ at percolation. Taking summation over all clusters of rank $r$ and mass $m$ we obtain: $$\begin{aligned}
\label{nrm0}
n_r(\rho_c)&=&\sum n_{r,m}(\rho_c)
=q_0\,\sum_{m_{\rm lo}}^{m_{\rm up}}\,m^{-\tau}\nonumber\\
&\sim&\frac{q_0}{\tau-1}
\left[(m_{\rm lo})^{-\tau+1}-(m_{\rm up})^{-\tau+1}\right]\nonumber\\
&=&\frac{q_0}{\tau-1}
\left[\left(\frac{m_{\rm lo}}{  m_r  }\right)^{-\tau+1}
-\left(\frac{m_{\rm up}}{  m_r  }\right)^{-\tau+1}\right]
m_r^{-\tau+1}.\end{aligned}$$ Our simulations suggest (not shown) that the mass distribution within a given rank is lognormal with a rank-independent standard deviation. Thus, for arbitrary upper and lower quantiles $m_{\rm up}$, $m_{\rm lo}$ of this distribution the values $$\frac{m_{\rm lo(up)}}{m_r}$$ are rank independent. Using this, we finally express $n_r(\rho_c)$ via $m_r$: $$\label{nrm}
n_r(\rho_c)=p_0\,  m_r  ^{-\tau+1}
\propto   m_r  ^{-1.05}.$$

The power law (\[nrm\]) is observed in a steady-state aggregation model of [@MNTG04] with index $1.147$. This index increase comparing to our $1.05$ is due to the fact that in [@MNTG04] intermediate clusters are removed from the lattice providing extra space for a larger number of smaller clusters.

Combining the mass-rank relation (\[MR\]) with (\[nrm\]) we obtain the following exponential rank distribution at percolation: $$\begin{aligned}
\label{GRr}
n_r(\rho_c)&\sim& p_0\,  m_r  ^{-\tau+1}
=p_0\,\left(c^{-\tau+1}\right)^{r-1}
= p_1\,10^{-b\,r}\end{aligned}$$ with $$p_1=p_0\,c^{\tau-1},~b=(\tau-1)\,\log_{10}c\approx 0.62.$$

This is indeed what we observe in Fig. \[fig\_GRr\] where the rank distribution $n_r$ at percolation is shown by the dash-dotted line. The study [@MNTG04] suggests $c^{1-\tau}=0.186$ while our predictions and observations lead to $c^{1-\tau}\approx 4.2^{-1.05} = 0.22$. The two values are in good agreement, the slight difference is explained, as in Eq. (\[nrm\]), by removal of intermediate clusters in [@MNTG04]. Next we consider the rank distribution for $\rho\ne\rho_c$.

Dynamical rank distribution: three-exponent scaling {#dynrankdist}
---------------------------------------------------

Here we expand results of the previous section by establishing the time-dependent rank dustribution. First, we consider the dynamics of rank population.

### Temporal dynamics of rank population

The dynamics of the total number $(n_r\cdot L^2)$ of the clusters of a given rank $r$ is illustrated in Fig. \[fig\_rank\] for $r=5,6,7$. The population follows a characteristic bell-shaped trajectory, with percolation at its rightward limb. As in the case of mass description, one does not observe steady-state behavior in the cluster dynamics: The population of each rank steadily develops to its peak as a result of merging of the clusters of lower ranks; then it starts decreasing, giving birth to the clusters of higher ranks. As naturally follows from the model definition, the peak of the population of a higher rank comes after the peak of a lower rank. Figure \[fig\_rank\_all\] shows the population dynamics for the ranks $1\le r \le 11$ in semilogarithmic scale. Here one clearly sees the similarity in the dynamics of different ranks. Note that this figure is remarkably similar to Fig. 7 from [@Tur+99] that shows the dynamics of clusters with logarithmically binned masses. We now proceed by establishing the appropriate time-dependent scaling laws.

### Time-dependent mass distribution

Recall that the temporal dynamics of the cluster mass distribution is given by the two-exponent scaling law [@Stauffer75; @SA; @Margolina+84]: $$\label{2exp}
n_m(\rho)\sim m^{-\tau}\,f_0(z),\qquad 
z=(\rho_c-\rho)m^{\sigma}+z_0,$$ with $\sigma=1/2$. The function $f_0$ has a bell-shaped form with maximum to the left of percolation; it can be roughly approximated by a Gaussian function [@Hoshen+79; @Margolina+84]: $$\label{f0}
f_0(z)\propto \exp\left(-a\,z^2\right).$$ Note that the shift $z_0$ is independent of $m$.

Considered as a function of $m$, the two-exponent scaling explains the power law mass distribution (\[GRm\]) at percolation (with $q_0=f_0(z_0)$) as well as the downward bend for $\rho<\rho_c$, clearly observed in Fig. \[fig\_GRm\] (dashed line); while as a function of $\rho$ it describes the bell-shaped dynamics of clusters with given mass $m$.

### Time-dependent rank distribution {#rscale}

Combining the scaling laws (\[MR\]) and (\[2exp\]) one formally obtains the two-exponent scaling for rank dynamics. However, the two exponent scaling does not work for ranks; to show this we assume more generally $$\label{r2exp}
n_r(\rho)\sim g_0(z)10^{-br},~
z=(\rho_c-\rho)h(r)+z'_0,$$ which is consistent with the exponential rank distribution of Eq.(\[GRr\]) at percolation (with $p_0=g_0(z'_0)$). Possible deviations from the pure exponential law at $\rho<\rho_c$ (clearly observed in Fig. \[fig\_GRr\]) and dynamics of a given rank (see Figs. \[fig\_rank\],\[fig\_rank\_all\]) are described by specific form of the functions $g_0(\cdot)$ and $h(\cdot)$. Following [@Hoshen+79] we define $$\label{nu}
\nu_r(z):=\frac{n_r(z)}{n_r(z'_0)}
=\frac{g_0(z)}{g_0(z'_0)}.$$ and choose $h(\cdot)$ in such a way that positions of the peaks of $\nu_r(z)$ coincide for different $r$; it is always possible by choosing the appropriate time change $h(r)$. Figure \[fig\_nu\]a shows the ratio $\nu_r(z)/\nu_1(z)$ for $r=2,3,6,8$. One can see that the two-exponent scaling does not work in our case: the curves do not coincide.

Nevertheless, the simple scaling picture is restored by introducing the additional, third, shift exponent: $$\label{r3exp}
h(r)= a_1\,10^{\sigma_1\,r},\qquad 
z'_0(r)=a_2\, 10^{-\sigma_2\,r}.$$ Function $g_0$ still can be approximated by a Gaussian function $$\label{g0}
g_0(z)\propto \exp\left(-\frac{z^2}{2}\right).$$

Once the correct scaling form is established, the use of (\[MR\]) is again legitimate, and the exponent $\sigma_1$ in Eq. (\[r3exp\]) can be evaluated as: $$\sigma_1=\sigma\log_{10} \hat{c}
\approx 0.24,$$ where $\hat{c}\approx 3$ is the median of $c$ values observed during $\rho<\rho_c$. The observed exponent $\sigma_1\approx 0.23$ (not shown) is fairly close to its predicted value. The shift exponent is estimated as $\sigma_2\approx 0.03$; while scale coefficients are $a_1\approx 1.54$, $a_2\approx 1.43$. The function $g_0(z)$ that uses these estimates is shown in Fig. \[fig\_g0\] where different symbols depict clusters of different ranks. The collapse is obvious, confirming the validity of the three-exponent scaling (\[r2exp\]), (\[r3exp\]), (\[g0\]).

In the scaling for cluster masses, the time renormalization $(\rho_c-\rho)m^{\sigma}$ collapses the dynamics of mass $m$ clusters onto the master curve $f_0(z-z_0)$ with its only peak shifted by $z_0$ leftward from percolation; the shift $z_0$ is mass independent. Similarly, in the scaling for ranks the time renormalization $(\rho_c-\rho)10^{\sigma_1\,r}$ collapses the dynamics of rank $r$ clusters onto the master curve $g_0(z-z'_0)$, although the shift now is rank dependent and is given by $10^{\sigma_2\,r}$. To illustrate this, we show the position of percolation on the righthand limb of the Gaussian $g_0(z-0.51)$ in Fig. \[fig\_nu\]b. The higher the rank, the closer the position of percolation to the peak of $g_0$.

Averaged scaling {#AS}
----------------

In applications, it is often impossible to measure the size distribution of system elements at a given time instant. Moreover, sometimes the instantaneous size distribution does not exist at all: This is indeed the case for the systems described by marked point processes widely used to model seismicity, volcano activity, starquakes, etc. [@DVJ]. In such situations one uses the averaged measurements. For instance, the famed Gutenberg-Richter law [@GR54; @Tur97; @BZ03] that gives exponential approximation to the size distribution of earthquakes (via their magnitudes) is valid only after appropriate averaging over a wide spatio-temporal domain. This explains the importance of the question: How do the distributions of Eq. (\[2exp\]), (\[r2exp\]) change after temporal averaging?

We answer this question for averaging over $0\le\rho\le\rho_c$. For the mass distribution this leads to: $$\begin{aligned}
\label{avem}
\widehat{n_m}&:=&\int_0^{\rho_c} n_m(\rho)d\rho
=\int_0^{\rho_c} f_0(z)\,m^{-\tau}d\rho\nonumber\\
&\propto&\int_0^{\rho_c}\exp\left\{-a\left[(\rho_c-\rho)m^{\sigma}
-z_0\right]^2\right\}m^{-\tau}d\rho\nonumber\\
&\propto& m^{-\tau-\sigma}\int_{u_1}^{u_2}
\exp\left\{u^2/2\right\}du\nonumber\\
&\propto& m^{-\tau-\sigma}~~(\approx m^{-5/2}).\end{aligned}$$ Here the last step neglects the weak dependence of the integral on $m$ (and uses the values $\tau\approx 2.0$, $\sigma=1/2$). The validity of (\[avem\]) is confirmed by the observed averaged mass distribution shown by the solid line in Fig. \[fig\_GRm\]. The averaged mass distribution is similar to that at percolation: it retains the power-law form while the slope is increased by $1/2$ due to averaging.

Similarly, we obtain for ranks: $$\begin{aligned}
\label{aver}
\widehat{n_r}&:=&\int_0^{\rho_c} n_r(\rho)d\rho
=\int_0^{\rho_c} g_0(z)\,10^{-br}d\rho\nonumber\\
&\propto&\int_0^{\rho_c}\exp\left\{-a'\left[
(\rho_c-\rho)10^{\sigma_1\,r}-a\,10^{-\sigma_2\,r}\right]^2
\right\}10^{-b\,r}d\rho\nonumber\\
&\propto& 10^{-(\sigma_1+b)\,r}\int_{u_1}^{u_2}
\exp\left\{u^2/2\right\}du\nonumber\\
&\propto& 10^{-(\sigma_1+b)\,r} = 
10^{r\,(1-\sigma-\tau)\log_{10}\tilde{c}}=
10^{-r\,\alpha_r}.\end{aligned}$$ The exponent $\alpha_r$ may vary from 0.71 to 0.93 depending on $3.0\le\tilde{c}\le 4.2$ (the range of $c$ values for the time when at least three ranks have been formed so the estimation of the distribution slope is meaningful). Simulations suggest (solid line in Fig. \[fig\_GRr\]) $\alpha_r=0.87$, which is in good agreement with our prediction. Again, the averaged rank distribution retains the exponential form of the distribution at percolation; while its index has increased due to averaging.

Correction to simple scaling {#corrections}
----------------------------

Due to finiteness of the lattice, the results of previous sections require some corrections to match exactly the simulated rank distributions. The appropriate corrections are described below.

### Corrected scaling at percolation

The pure power and exponential laws in Figs. \[fig\_GRm\], \[fig\_GRr\] are just first-order approximations to the observed cluster distributions at percolation. In both cases one sees the downward bending for small clusters and upward bending for large clusters. These are not due to statistical fluctuations. The downward bending for small clusters is explained by “deviations from scaling” [@Hoshen+79]: it can be shown analytically that the small clusters do not yet obey the general scaling law of Eqs. (\[GRm\]), (\[GRr\]) which holds only for large enough masses (ranks). The upward bend at large clusters is due to finite-size effects [@Hoshen+79; @Margolina+84]: each large cluster that reaches outside the lattice boundary is “seen” as a number of smaller clusters, thus creating the upward deviation from the pure power (exponential) law. This phenomenon is especially important when the system is close to percolation and clusters of arbitrary large sizes have already been formed. The appropriate scale corrections for the mass distribution were studied by Hoshen et al. [@Hoshen+79] and Margolina et al. [@Margolina+84].

To study the above phenomena it is convenient to consider the normalized functions $$N_m:=m^{\tau-1}\sum_{m'\ge m} n_{m'},\qquad
N_r:=10^{br} n_r,$$ which, in the absence of scale corrections, would become constants: $$N_m=\frac{q_0}{\tau-1},\qquad
N_r=p_1\,c^{\tau-1}.$$ The function $N_r$ is shown in Fig. \[fig\_Nr\]a; it clearly deviates from the horizontal plateau at both sides.

In case of the mass distribution, the corrections to scaling are given by [@Margolina+84]: $$\label{GRmc}
n_m(\rho_c)\simeq m^{-\tau}
\left(q_0+q_1\,m^{-\Omega}+q_L\,m^{1/D}L^{-1}\right),$$ where $\Omega \approx 0.75$, $1/D=48/91$ is the universal mean cluster radius exponent, and $q_0,q_1,q_L$ are independent of $s$ and $L$. The first additional term describes the deviation from scaling for small clusters, while the second one is responsible for finite-size effects.

For rank distribution, the “deviations from scaling” at lower clusters are only observed for $r=1$; while the finite-size effects at large clusters are clearly present for many ranks. Accordingly, we propose the following correction to scaling for the rank distribution: $$\label{GRrc}
n_r(\rho_c)\simeq 10^{-br}
\left(p_0+p_L\,10^{d\,r}L^{-1}\right),~r>1.$$ with $$d=\frac{1}{D}\,\log_{10}c \approx 0.33.$$ The observed value of $d$ can be estimated by plotting $(n_r\,10^{br}-p_0)$ as a function of $r$ as shown in Fig. \[fig\_Nr\]b. The observed ranks $4\le r \le 9$ follow the predicted scaling (\[GRrc\]) nicely.

Importantly, the corrections to scaling (\[GRrc\]) act at all cluster sizes, so they can not be neglected even for the intermediate clusters, not only for the largest ones. Indeed, their effect decreases with $L$, but this decrease is very slow. Notably, as shown by Morein et al. [@MNTG04] (their Fig. 5) even for lattices as large as $L=30,000$ during the process when clusters as large as 2% of the lattice size are removed, the cluster size distribution clearly exhibits the upward deviations at large ranks ($r=11,12,13$.) For smaller systems these deviations become dominant and may lead to an artificial decrease of the observed slope of cluster size distribution; this is demonstrated in Fig. \[fig\_GRm\],\[fig\_GRr\] and is also seen in the analysis of Turcotte et al. [@Tur+99] (their Fig. 9).

### Dynamics of scaling corrections {#scalecorr}

Since the finite size effects play an important role in shaping the observed cluster size distribution, it is worth studying their dynamics. Specifically, we will be interested in transition of the cluster size distribution from the convex shape (in semi- or bilogarithmic scale) at $\rho\ll\rho_c$ to formation of the upward bend at percolation.

For this we introduce a measure of convexity for the rank distribution, defined as an area between $\log n_r(\rho)$ and a chord connecting its first and last points as shown in Fig. \[fig\_mu\] (the point $r=1$ is not considered being affected by the deviations from scaling): $$\label{mconv}
\mu(\rho):=\int_{2}^{r_{\rm max}} 
\left[\log_{10} n_r(\rho)-(A\,r+B)\right]
\,d\,r,$$ with $$A=\frac{\log_{10}\left(n_{r_{\rm max}}/n_2\right)}
{r_{\rm max}-2},
~B=\log_{10}\,n_2-2A.$$ The values of $\mu$ are positive when $n_r(\rho)$ is convex in semilogarithmic scale, negative when it is concave, and vanish when it is linear. The measure $\mu(\rho)$ averaged over 1,000 runs on the lattice $L=2000$ is shown in Fig. \[fig\_mu\]; the bell-shaped form of $\mu$ is decorated by the logperiodic oscillations for $(\rho_c-\rho)>10^{-2}$ explained by creation of new ranks, which temporarily increases convexity. Zero level is crossed at about $(\rho_c-\rho)=2\cdot 10^{-3}$, after that the rank distribution is concave. A detailed analysis (not shown) demonstrates that the distribution is never exactly linear; the transition from convex to concave shape is realized through the wave-shaped form when the distribution is still convex for the lower $r$, but is already concave for the higher ones. Qualitatively the same picture is observed for the mass distribution $n_m(\rho)$ (in bilogarithmic scale).

The transformation of the cluster size distribution prior to percolation is not unlike a well-known pattern “upward bend” first described by Narkunskaya and Shnirman [@NS90; @NS94] in an early static model of defect development. Later it was found in steel samples and seismicity of California [@RKB97], and confirmed by the dynamical modeling of failure in a hierarchical system (so-called colliding cascade models) [@GKZN00; @ZKG03].

Mass dynamics of a given rank {#MDGR}
-----------------------------

Here we consider the dynamics of total and average mass of rank $r$ clusters: $$\label{mass}
M_r=\sum m\,n_{rm},\qquad
  m_r   =\frac{\sum m\,n_{rm}}{\sum n_r}
=\frac{M_r}{n_r}.$$ Here $n_{rm}$ denotes the number of clusters of rank $r$ and mass $m$. Figure \[fig\_avem\] shows $n_r$, $M_r$, and $  m_r   $ for rank 5; the similar picture is observed for other ranks. It is tempting to use Gaussian approximation for $M_r$ and predict the Gaussian dynamics of $  m_r   $ (as a ratio of two Gaussians) and relate their parameters. Detailed analysis however demonstrates that under this approach the peak of $  m_r   $ for ranks $r\ge 9$ should be observed after percolation; while in simulations this peak is always prior to percolation (not shown). Note that one still might approximate $M_r$ and $  m_r   $ by Gaussians with properly scaled parameters; such approximations will be good for rough curve fitting, but will fail to reproduce deeper properties of cluster dynamics. This demonstrates the general limitations of Gaussian approximations in the percolation problem.

Cluster fractal structure {#CFS}
=========================

In this section we evaluate the fractal structure of clusters considering the mass-circumference relation $$\label{circum}
m \propto C^{D_r},$$ where $C$ is the number of empty neighbors of a cluster of mass $m$. For percolation cluster at infinite grid we have $D_r=1$, which shows that the percolation cluster is a “linear” rather than a space-filling object [@SA]. Figure \[fig\_D\]a shows the cluster masses as a function of their circumference for different ranks. It is easily seen how the linear scaling $D_r=1$ gradually develops as rank increases. Figure \[fig\_D\]b shows the index $D_r$ estimated for $1\le r\le 9$.

Figure \[fig\_fin\] shows the dynamics of $D_5$ prior to percolation; noteworthy, its steady state behavior is altered by a gradual increase as $\rho\to\rho_c$. Similar increase is observed for clusters of other ranks.

To explain the increase of $D_r$ we recall that the rate of cluster coalescence is proportional to their circumference (see e.g. [@GNT99]). Thus, for a given mass, clusters with a lower $D_r$ have larger circumference, and a higher chance to coalesce. When a sufficient number of rank $r$ clusters have been formed, the clusters with low $D_r$ start to coalesce leaving the high-$D_r$ clusters on the grid.

Another reason for the increase of $D_r$ is the finite-size effects. Specifically, this is an effect of having clusters that on an infinite grid have already gained higher ranks, but on our finite lattice are still small.

Dynamical constraint {#DC}
====================

Here we report an interesting regularity in rank dynamics that put a notable constraint on analytical modeling of percolation process. Specificaly, we consider the slope between two consecutive points of the rank distribution: $$\theta_r(\rho):=\log\,\frac{n_r(\rho)}{n_{r+1}(\rho)}.$$ Dynamics of $\theta_4$ is shown in Fig. \[fig\_peaks\]a together with that of $n_6$. Noteworthy, the peaks of two curves (minimum of $\theta_4$ and maximum of $n_6$) coincide. This happens to be true for all ranks: positions of corresponding peaks are shown as a function of rank in Fig. \[fig\_peaks\]b. Such perfect matching is very unlikely to be accidental. Thus we conjecture that in order for $n_r(\rho)$ to properly describe the time-dependent behavior of rank population, the following system of differential equations must have a solution: $$\label{constraint}
\left\{\begin{array}{l}
\dot{\theta}_r=0\\
\dot{n}_{r+2}=0
\end{array}
\right.
\Rightarrow
\left\{\begin{array}{l}
\dot{n_r}\,n_{r+1}-n_r\dot{n}_{r+1}=0\\
\dot{n}_{r+2}=0
\end{array}
\right.$$

Applying this constraint to the three-exponent scaling of Eqs. (\[r2exp\]),(\[r3exp\]),(\[g0\]) we find $$\label{sigmas}
\sigma_2=\sigma_1+\log_{10}\left(1-10^{-2\sigma_1}\right).$$ According to (\[sigmas\]), the observed value $\sigma_1=0.23$ gives $\sigma_2=0.04$, which is 33% larger than the observed value $\sigma_2=0.03$. The discrepancy is due to the approximate character of the Gaussian approximation (\[g0\]) for $g_0$.

Discussion
==========

The goal of this study was to describe the evolution of percolation model in terms of consecutive aggregation of smaller clusters into larger ones using the Horton-Strahler hierarchical scheme. First, this contributes to a novel understanding of the percolation phenomenon as a time-dependent hierarchical inverse cascade process. Second, this allows one to test the validity of the approach introduced by Gabrielov et al. [@GNT99] and further developed by Morein et al. [@MNTG04] for a steady-state approximation to a general aggregation process.

We considered dynamical and scaling properties of site-percolation on a 2[*D*]{} square lattice. Following [@GNT99] we described clusters by hierarchical trees that reflect the history of cluster formation; the Horton-Strahler scheme was used to rank the trees and thus the corresponding clusters. We concentrated on the development of the first percolation cluster, thus working with a system that does not exhibit the steady-state dynamics, contrary to the studies [@GNT99; @MNTG04] that have developed mean-field steady-state approximations to the system.

Combining the results obtained in the classical percolation studies with the Tokunaga constraint on the cluster branching structure we derived various rank-dependent scaling laws connecting the number $n_r$ of clusters of rank $r$, their average mass $m_r$, and the rank $r$. We have compared the parameters of these laws with those predicted and observed by [@GNT99; @MNTG04] in steady-state aggregation models. The values of parameters are shown to be in a perfect agreement, confirming the validity of the approach used in [@GNT99; @MNTG04]. In absence of the steady-state behavior, we derived the time-dependent versions of the scaling laws. We reported the three-index scaling (\[r2exp\]), (\[r3exp\]) for the number $n_r(\rho)$ of clusters of a given rank, which deviates from the classical two-exponent scaling for masses.

We studied in detail the transition of the system from earlier stages to the vicinity of percolation and reported several characteristic phenomena observed as $\rho\to\rho_c$. They include transformation of the cluster size distribution not unlike that observed in seismicity, steel samples, and previous models of hierarchical fractures [@NS90; @RKB97; @GKZN00; @ZKG03]; and increase of the cluster fractal dimension. In our simple model these phenomena are partly explained (qualitatively as well as quantitatively) by finite-size effects; nevertheless we believe that they should not be neglected as irrelevant side-effects of numerical simulation. In fact, in practice we often work with systems that are described by intermediate depth hierarchies (in other words they have intermediate number of degrees of freedom). The percolation results related to small and intermediate lattices might be of high relevance in describing such systems. In addition, simulations on large lattices ($L=30,000$) performed by Morein et al. [@MNTG04] show that finite-size effects are still present even for large $L$.

We have formulated an empirical constraint of Eq. (\[constraint\]) for the time-dependent behavior of rank population size $n_r(\rho)$; the constraint follows very clearly from the observed values of $n_r(\rho)$. It would be interesting to check this condition in real systems traditionally described by the percolation model.

Our closing remark is on the index $\tau$ of cluster mass distribution at percolation (Eq. (\[GRm\])). The studies of Gabrielov et al. [@GNT99] and Morein et al. [@MNTG04] predict $\tau=2$; which slightly deviates from the well established theoretical value of the Fisher exponent $\tau=187/91\approx 2.05$. The index of the mass distribution is an essential characteristic of a system, thus even this slight difference of 2.5% might seem disappointing implying the intrinsically approximate character of the modeling of [@GNT99; @MNTG04]. In fact, this implication is not true. To validate the approach of [@GNT99; @MNTG04] we notice that the Fisher exponent is tightly connected to the precise count of cluster particles on a site-level, hardly feasible in practice. At the same time, the studies [@ZLK99; @CZ03] have demonstrated that when we “characterize the size distribution of clusters in a way that circumvents the site-level description” considering any “macroscopic measure of the length scale of the cluster”, the exponent of the corresponding scaling law becomes $2$, universally for all 2[*D*]{} systems. An example of a “macroscopic measure” is the linear size in arbitrary direction, the radius of gyration, the diameter of the covering disk, etc. Clearly, the modeling of [@GNT99; @MNTG04] deals with such a macroscopic measure of cluster size, and hence predicts the correct scaling exponent.

We are sincerely grateful to Bill Newman for numerous focused discussions; his critics and advice have helped significantly in organizing and presenting the results. We thank Vladimir Keilis-Borok, Gleb Morein, and Donald Turcotte for their continuous interest to the work and David Shatto for help in preparing the manuscript. This study was partly supported by NSF, grant ATM 0327558.

[99]{}

Allegre, C.J., LeMouel, J.L., Provost, A. (1982) Scaling rules in rock fracture and possible implications for earthquake prediction. [*Nature*]{}, [**297**]{}, (5861): 47-49.

Badii, R., and Politi, A. (1997) [*Complexity: Hierarchical Structures and Scaling in Physics.*]{} Cambridge University Press. 318 pp.

Ben-Zion, Y. (2003) Appendix 2, Key Formulas in Earthquake Seismology, in [*International Handbook of Earthquake and Engineering Seismology*]{}, Part B, 1857-1875, Academic Press.

Blanter, E.M., Shnirman, M.G., LeMouel, J.L., Allegre, C.J. (1997) Scaling laws in blocks dynamics and dynamic self-organized criticality. [*Phys. Earth Planet. Inter.*]{} [**99**]{} (3-4): 295-307.

Blanter, E.M., Shnirman, M.G., LeMouel, J.L. (1997) Hierarchical model of seismicity: scaling and predictability. [*Phys. Earth Planet. Inter.*]{}, [**103**]{} (1-2): 135-150.

Cardy, R., Ziff, R.M. (2003) Exact Results for the Universal Area Distribution of Clusters in Percolation,Ising,and Potts Models. [*J. Stat. Phys.*]{}, [**110**]{} (1-2): 1-33.

da Costa, F.P., Grinfeld, M., Wattis, J.A.D. (2002) A hierarchical cluster system based on Horton-Strahler rules for river networks. [*Studies Appl. Math.*]{} [**109**]{} (3): 163-204.

Daley, D.J., Vere-Jones, D. (2003) [*An Introduction to the Theory of Point Processes. Vol. I,*]{} 2nd ed., Springer, New-York, pp. 469.

Fisher, M. E. (1967) [*Physics*]{} [**3**]{}, 255.

Gabrielov, A., Keilis-Borok, V., Zaliapin, I., Newman, W.I. (2000) Critical transitions in colliding cascades. [*Phys. Rev. E*]{}, [**62**]{}, 237-249.

Gabrielov, A., Newman, W.I., Turcotte, D.L. (1999) An exactly soluble hierarchical clustering model: inverse cascades, self-similarity, and scaling. [*Phys. Rev. E*]{}, [**60**]{}, 5293-5300.

Gutenberg, B., Richter, C.F. (1954) [*Seismicity of the Earth and Associated Phenomena*]{}, Princeton University Press, Princeton, 2-nd ed.

Horton, R.E., (1945) Erosional development of streams and their drainage basins: Hydrophysical approach to quantitative morphology, [*Geol. Soc. Am. Bull.*]{}, [**56**]{}, 275-370.

Hoshen, J., Stauffer, D., Bishop, G.H., Harrison, R.J., and Quinn, G.D. (1979) Monte Carlo experiments on cluster size distribution in percolation. [*J. Phys. A*]{}, [**12**]{}, 1285-1303.

Knopoff, L., Newman, W.I. (1983) Crack fusion as a model for repetitive seismicity. [*Pure Appl. Geophys.*]{}, [**121**]{}, (3): 495-510.

Leyvraz, F. (2003) Scaling theory and exactly solved models in the kinetics of irreversible aggregation. [*Phys. Rep.*]{}, [**383**]{}, (2-3): 95-212.

Margolina, A., Nakanishi, H., Stauffer, D., and Stanley, H.E. (1984) Monte Carlo and series study of corrections to scaling in two-dimensional percolation. [*J. Phys. A*]{}, [**17**]{}, 1683-1701.

Morein, G., Newman, W.I., Turcotte, D.L., and Gabrielov A. (2004) An inverse cascade model for self-organized complexity and natural hazards. Preprint, pp. 9.

Narkunskaya, G. S. and Shnirman, M. G., (1990) Hierarchical model of defect development and seismicity. Phys. Earth. Planet. Inter., [**61**]{}, 29-35.

Narkunskaya, G. S. and Shnirman, M. G., (1994) An algorithm of earthquake prediction. Computational Seismology and Geodynamics (AGU, Washington, D.C.), [**1**]{}, 20-24.

Newman, W.I., Gabrielov, A. (1991) Failure of hierarchical distributions of fiber bundles. I. [*Internat. J. of Fracture*]{}, [**50**]{}, 1-14.

Newman, W.I., Knopoff, L. (1982) Crack fusion dynamics – a model for large earthquakes. [*Geophys. Res. Lett.*]{}, [**9**]{} (7): 735-738.

Newman, W.I., Knopoff, L. (1983) A model for repetitive cycles of large earthquakes. [*Geophys. Res. Lett.*]{}, [**10**]{} (4): 305-308.

Newman, W.I., Knopoff, L. (1990) Scale-invariance in brittle-fracture and the dynamics of crack fusion. [*Intl. J. of Fracture*]{}, [**43**]{} (1): 19-24.

Newman, M.E.J., Ziff, R.M. (2001) Fast Monte Carlo algortithm for site or bond percolation. [*Phys. Rev. E*]{}, [**64**]{}, 016706.

Newman, W.I., Turcotte, D.L., and Gabrielov, A. (1997) Fractal trees with side branching. [*Fractals*]{}, [**5**]{}, 603-614.

Newman, W.I., Turcotte, D.L., and Gabrielov, A. (1995) Log-periodic behavior of a hierarchical failure model with applications to precursory seismic activation. [*Phys. Rev. E*]{}, [**52**]{}, 4827-4835.

Rotwain, I., Keilis-Borok, V. and Botvina, L. (1997) Premonitory transformation of steel fracturing and seismicity, [*Phys. Earth Planet. Inter.*]{}, [**101**]{}, 61-71.

Shnirman, M.G., Blanter, E.M. (2001) Criticality in a dynamic mixed system. [*Phys. Rev. E*]{}, [**64**]{} (5): Art. No. 056123 Part 2.

Sornette, D., Johansen, A., Arneodo, A., Muzy, J.-F., and Sauleur, H. (1996) Complex fractal dimensions describe the internal hierarchical structure of DLA, [*Phys. Rev. Lett.*]{} [**76**]{}, 251-254.

Sornette, D. (2004) [*Critical Phenomena in Natural Sciences.*]{} 2-nd ed., Springer-Verlag Heidelberg, pp. 528.

Stauffer, D. (1975) Violation of dynamical scaling for randomly dilute Ising ferromagnets near percolation threshold. [*Phys. Rev. Lett.*]{} [**35**]{} (6): 394-397.

Stauffer, D., Sornette, D. (1998) Log-periodic oscillations for biased diffusion on random lattice, [*Physica A*]{}, [**252**]{} (3-4): 271-277.

Stauffer, D., Aharony, A. (1994) [*Introduction to percolation theory.*]{} 2-nd ed., Taylor & Francis.

Strahler A. N. (1957) Quantitative analysis of watershed geomorphology. [*Trans. Am. Geophys. Un.*]{}, [**38**]{}, 913-920.

Tokunaga, E. (1978) Consideration on the composition of drainage networks and their evolution. [*Geographical Rep. Tokyo Metro. Univ.*]{}, [**13**]{}, 1-27.

Toroczkai, Z. (2001) Topological classification of the Horton-Strahler index on binary trees, [*Phys. Rev. E*]{}, [**65**]{}, 016130.

Turcotte, D.L. (1997) [*Fractals and Chaos in Geology and Geophysics.*]{} (2nd edition). Cambridge University Press. 398 pp.

Turcotte, D.L., Malamud, B.D., Morein, G., and Newman W.I. (1999) An inverse cascade model for self-organized critical behavior, [*Physica A*]{}, [**268**]{}, 629-643.

Turcotte, D.L., Malamud, B.D., Guzzetti, F., Reichenbach, P. (2002) Self-organization, the cascade model, and natural hazards. [*Proc. Natl. Ac. Sci.*]{}, [**99**]{}, 2530-2537.

Turcotte, D.L., Pelletier, J.D., Newman, W.I. (1998) Networks with side branching in biology. [*J. Theor. Biology*]{}, [**193**]{} (4): 577-592.

Ziff, R.M., Lorenz, C.D., Kleban, P. (1999) Shape-dependent universality in percolation [*Physica A*]{}, [**266**]{} (1-4): 17-26.

Zaliapin, I., Keilis-Borok, V., Ghil, M. (2003) A Boolean delay equation model of colliding cascades. Part II: Prediction of critical transitions. [*J. Stat. Phys.*]{} [**111**]{}, (3-4): 839-861.

![Multiple coalescence of clusters. a) Coalescence of clusters is materialized by adding to the lattice a new particle [**N**]{} (black) that is a neighbor to one, two, three, or four existing clusters (numbered gray sites). The relative frequencies $Q_k$ , $k=1,2,3,4$ of $k$-coalescences based on similations with $L=2,000$ are shown in the figure. The corresponding tree is constructed as shown in panel b) (for $k=1$) and c) (for $k=2$). The cases $k=3,4$ are analogous to $k=2$. Note that about 95% of coalescences result in merging two clusters. See text for details.[]{data-label="fig_model1"}](figure1.eps)

![Tree representation of clusters: scheme. The dynamics is from left to right. At first step particle [**A**]{} is dropped onto the lattice and a one-particle cluster is formed; it is represented by a one-node tree. At second step another one-particle cluster [**B**]{} is formed; it is represented by another one-node tree. At third step new particle [**C**]{} coalesces with cluster [**A**]{} to form two-particle cluster [**AC**]{}. This cluster is represented by a three-node tree; note that the connecting node of the tree does not correspond to any particle. At fourth step new particle [**D**]{} connects existing clusters [**AC**]{} and [**B**]{} to form four-particle cluster [**ABCD**]{}. This cluster correspond to a five-node tree.[]{data-label="fig_model2"}](figure2.eps)

![a) Non-uniqueness of tree representation. Two different 11-particle clusters that correspond to the same tree shown in panel b). Particles have been dropped according to their alphabet marks; so first was the particle [**A**]{}, then [**B**]{}, etc. b) Horton-Strahler ranking: illustration. The ranks are shown next to the tree nodes.[]{data-label="fig_model3"}](figure3.eps)

![Mass-rank distribution observed on a 2,000$\times$2,000 lattice at percolation. Dots – individual clusters, balls – average mass $m_r$ within a given rank. Line shows the relation $m_r=\left[10^{\,0.625}\right]^{\,r-1}=4.2^{\,r-1}$.[]{data-label="fig_MR"}](figure4.eps)

![Parameter $c$ of the mass-rank relation $m_r=c^{r-1}$ as a function of time. At percolation $c(\rho_c)\approx 4.2$; the Euclidean limit of [@GNT99] corresponds to $c=3.25$, it is reached at $\rho_c-\rho\approx 0.14$.[]{data-label="fig_c"}](figure5.eps)

![Mass distribution of clusters observed on a $2,000\times 2,000$ lattice at percolation $\rho=\rho_c$ (dash-dotted line), $\rho=0.48$ (dashed line), and averaged over $0<\rho<\rho_c$ (solid line). To smooth out statistical fluctuations we show the cumulative distribution: $\propto\sum_{m'>m}n_m$. For comparison, all curves are normalized to unity at $m=1$.[]{data-label="fig_GRm"}](figure6.eps)

![Rank distribution of clusters observed for $2,000\times 2,000$ lattice at percolation $\rho=\rho_c$ (dash-dotted line), $\rho=0.29$ (dashed line), and averaged over the percolation cycle $0<\rho<\rho_c$ (solid line). For comparison, all curves are normalized to unity at $r=1$.[]{data-label="fig_GRr"}](figure7.eps)

![Dynamics of population $n_r\cdot L^2$ of a given rank, $r=5,6,7$ for $L=2,000$. Moment of percolation is depicted by a vertical dashed line.[]{data-label="fig_rank"}](figure8.eps)

![Dynamics of population $n_r$ of a given rank, $1\le r\le11$ in semilogarithmic scale. Moment of percolation is depicted by a vertical dashed line. (Cf. Fig. 7 from [@Tur+99]).[]{data-label="fig_rank_all"}](figure9.eps)

![Scaling for rank dynamics. a) Ratios $\nu_r(z)/\nu_1(z)$ do not collapse thus rejecting the two-exponent scaling hypothesis; see details in Sect. \[rscale\]. b) Position of percolation on the normalized Gaussian $g_0(z-0.51)$; see details in Sect. \[rscale\].[]{data-label="fig_nu"}](figure10.eps)

![Three-exponent scaling for rank dynamics. The master Gaussians $g_0(z)$ for different ranks collapse when using the renormalization given by Eqs. (\[r2exp\]),(\[r3exp\]).[]{data-label="fig_g0"}](figure11.eps)

![Corrections to scaling. The pure exponential rank distribution of Eq. (\[GRr\]) suggests a horizontal plateau for the normalized function $N_r=10^{br}n_r$, while the observed values clearly deviate from the plateau at small and large clusters (panel a). The large cluster deviation is due to finite size efffects and is described by an exponential correction of Eq. (\[GRrc\]) with $d\approx 0.33$ (panel b).[]{data-label="fig_Nr"}](figure12.eps)

![Dynamics of scale corrections. A convexity measure $\mu(\rho)$ is defined by Eq. (\[mconv\]) and illustrated in the insert. It is positive for convex, and negative for concave rank distribution. The downward bend of the rank distribution observed at early stages ($\mu>0$) is changed to the upward one ($\mu<0$) for $(\rho_c-\rho)<2\cdot10^{-3}$. See details in Sect. \[scalecorr\].[]{data-label="fig_mu"}](figure13.eps)

![Dynamics of number of clusters $n_r$ (solid line), total mass $M_r$ (dashed line), and average mass $m_r$ (dash-dotted line) for clusters of rank $r=5$.[]{data-label="fig_avem"}](figure14.eps)

![Fractal structure of clusters. a) Mass-circumference relation for clusters of different ranks. The asymptotic power relation with slope 1.0 is gradually develops as rank increases. b) Values of fractal dimension $D_r$ (Eq. (\[circum\])) for different ranks. Both panels correspond to a $2,000\times2,000$ lattice at percolation.[]{data-label="fig_D"}](figure15.eps)

![Premonitory increase of cluster fractal dimension. The steady-state dynamics of fractal dimension $D_5$ (Eq. \[circum\]) changes, and $D_5$ starts to increase, as system approaches percolation. Similar phenomenon is observed for other ranks. []{data-label="fig_fin"}](figure16.eps)

![Dynamical constraint for $n_r(\rho)$. Dynamics of $\theta_4=\log(n_4/n_5)$ and $n_6$ is shown in panel a: peaks of two curves coincide. The similar phenomenon is observed for other ranks: panel b shows the times of maxima of $n_r$ (circles) and minima of $\theta_{r-2}$ (triangles) for $3\le r\le 9$.[]{data-label="fig_peaks"}](figure17.eps)
