Hierarchical aggregation in percolation model

Tectonophysics 413 (2006) 93 – 107 www.elsevier.com/locate/tecto

Hierarchical aggregation in percolation model Ilya Zaliapin a,*, Henry Wong a, Andrei Gabrielov b a

Institute of Geophysics and Planetary Physics, University of California Los Angeles, Slichter Hall 3845, CA, 90095-1567, USA b Departments of Mathematics and Earth and Atmospheric Sciences, Purdue University, West Lafayette, IN, 47907, USA Received 27 December 2004; received in revised form 18 May 2005; accepted 4 October 2005 Available online 28 November 2005

Abstract There is a growing belief that the complex dynamics of seismicity can be better understood by studying the collective behavior of numerous lithosphere instability sources rather than focusing on the details of each of them. Classical site-percolation is a simple and tractable model which exhibits such important general features of complex systems as criticality and phase transitions of second kind. It also illustrates the mechanism of hierarchical aggregation, which is very important for explaining collective phenomena in material fracture and earthquake nucleation processes. We study the dynamics of a 2D site percolation model on a square lattice using the hierarchical approach introduced by Gabrielov et al., Phys. Rev. E., 5293–5300, 1999. The key elements of the approach are the tree representation of clusters and the Horton–Strahler scheme for cluster ranking. Accordingly, the evolution of percolation model is considered as a hierarchical inverse cascade of cluster aggregation. We analyzed the growth of the percolation cluster and established the time-dependent rank distribution of its subclusters, as well as corresponding laws for its mass, rank, and their relationship. We report several phenomena premonitory to the onset of percolation that complement the traditional power-law increase of the model’s observables. In addition, we have shown that the Tokunaga side-branching constraint uniquely determines the mass–rank relationship for a general aggregation process (not necessarily originated from the percolation model). The results can be used for development and improvement of earthquake prediction techniques. D 2005 Elsevier B.V. All rights reserved. Keywords: Hierarchical aggregation; Percolation; Earthquake prediction

1. Introduction Earthquakes pose an intolerable threat to society; yet their complex dynamics can hardly be described by a constitutive set of a few fundamental equations. There is an increasing belief that predictive understanding of the behavior of the lithosphere, which generates earthquakes, might result if the research focus shifts from an increasing number of specific sources of lithospheric instability and deformation to their essential collective * Corresponding author. E-mail addresses: [email protected] (I. Zaliapin), [email protected] (H. Wong), [email protected] (A. Gabrielov). 0040-1951/$ - see front matter D 2005 Elsevier B.V. All rights reserved. doi:10.1016/j.tecto.2005.10.010

behavior (Aki, 1995; Keilis-Borok, 1996, 2002; KeilisBorok and Soloviev, 2003; Turcotte, 1991, 1997, 2001; Rundle et al., 2000). During the last two decades, several general concepts and approaches for studying collective phenomena have been successfully applied to various geophysical problems, particularly to describing the dynamics of seismicity. Among them are the concepts of phase transition of second kind (Kadanoff, 2000; Rundle et al., 2000) and self-organized criticality (SOC) (Bak et al., 1998). The latter applies to systems that are maintained near a critical point: a state with no natural length scale typically reflected in a power-law statistics of a system’s observables (think of the powerlaw distribution of seismic moments).

94

I. Zaliapin et al. / Tectonophysics 413 (2006) 93–107

Classical site-percolation model presents probably the simplest and best studied system that exhibits criticality and experiences (geometrical) phase transition of second kind (Stauffer and Aharony, 1994). It is widely used as a toy model for spatially distributed stochastic processes, such as diffusion in disordered media, forest fires, gelation, semiconduction, etc. Importantly for our study, percolation model presents a very transparent mechanism of the process of hierarchical aggregation, which is very relevant for the dynamics of seismicity. An early idea of using hierarchical aggregation for describing essential properties of material fracture and earthquake nucleation appeared in the seminal works of Allegre et al. (1982), and Newman and Knopoff (1982, 1983, 1990), Knopoff and Newman (1983) and has been actively employed since then (Blanter et al., 1997a,b; Gabrielov et al., 2000a; Narkunskaya and Shnirman, 1990; Newman and Gabrielov, 1991; Newman et al., 1995; Shnirman and Blanter, 2001; Zaliapin et al., 2003). A recent general review of the theory and models of kinetics of irreversible aggregation is given by Leyvraz (2003). The principal novelty facilitating the recent efforts in the hierarchical aggregation research is the use of the Horton–Strahler scheme (Horton, 1945; Strahler, 1957; Newman et al., 1997) to assign ranks to the elements of hierarchical structures (cascades), as well as the Tokunaga classification (Tokunaga, 1978; Newman et al., 1997), which is used to describe the complex selfsimilar branching within those structures. The rediscovery of the Horton–Strahler scheme (see Section 3), which originated in hydrology and has not been well known in physical applications, allowed Gabrielov et al. (1999) to formulate an exactly solvable model for steady-state approximation to the aggregation process. Their analytical results were confirmed by extensive simulations of Morein et al. (2004), who focused on a steady-state regime in a modified forest-fire model; and Zaliapin et al. (2004), who applied the hierarchical description to the classical percolation model. The Horton–Strahler ranking scheme is shown essential for building an empirical description and developing numerical models in various geophysical, biological, and computational applications (da Costa et al., 2002; Badii and Politi, 1997; Gabrielov et al., 1999; Morein et al., 2004; Newman et al., 1997; Toroczkai, 2001; Turcotte et al., 1998). Notably, hierarchical ranks are more feasible than masses for observations in practice; they present a convenient bmacroscopicQ measure of cluster size, and provide a natural alternative to the logarithmic binning commonly used in analysis of modeled and observed processes.

In this study we describe the evolution of a percolation model (and, particularly, the percolation cluster) in terms of hierarchical aggregation of smaller clusters into larger ones using the Horton–Strahler ranking scheme. Our goal is three-fold: First, this contributes to a novel understanding of the percolation phenomenon as a time-dependent hierarchical inverse cascade process (Turcotte et al., 1999). Second, this further validates the analytical modeling approach introduced by Gabrielov et al. (1999) and developed by Morein et al. (2004) for a steady-state approximation to a general aggregation process. Third, this sets a basis for developing earthquake prediction techniques based on dynamical properties of the inverse cascades. A simple model considered in this work naturally suggests several patterns premonitory to the percolation onset; they could be further tested in this and other models and observations. Specifically, we describe the evolution of percolation model in terms of an inverse cascade of hierarchical cluster aggregation. We focus on the evolution of the first spanning cluster in the classical site-percolation model, and describe it by tracing the 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 (Gabrielov et al., 1999; Morein et al., 2004; Turcotte et al., 1999) as well as from the asymptotic focus on the percolation onset typical for the classical percolation studies (Stauffer and Aharony, 1994). Many phenomena encountered in the percolation model mimic what we see when the phase transitions of second kind occur. Note however that these phenomena are of purely geometrical and statistical rather than physical nature. Indeed, the physical percolation theory is largely 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. 2. Model 2.1. Dynamics We consider the classical 2D site-percolation model (Stauffer and Aharony, 1994). The model dynamics

I. Zaliapin et al. / Tectonophysics 413 (2006) 93–107

Fig. 1. Sample percolation model. The percolation cluster is shown in black: it connects the top and bottom sides of the grid. Coincident smaller clusters are shown in gray, empty sites are white.

starts with an empty L
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 (Stauffer and Aharony, 1994). 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 q = 0 at the start to q = 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 (Fig. 1). In the limit of an infinite lattice, the percolation cluster is formed at q c c 0.59274606 (Newman and Ziff, 2001), while for a finite lattice the critical time q c (L) is smaller (Stauffer and Aharony, 1994): qc ð LÞ ¼ qc  cL3 :

95

particle which will be a neighbor to one or more existing clusters. Fig. 2a illustrates the four possible types of coalescence; we call k-coalescence a situation when a newly dropped particle (marked N in the figure) is a neighbor to exactly k = 1, 2, 3 or 4 existing clusters (gray numbered sites). Numerical simulations on a square lattice with L = 2000 suggest the following relative frequencies Q k of k-coalescences: Q 1 c 0.628, Q 2 c0.318, Q 3 c 0.052, Q 4 c 0.002. Fig. 2b,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 a neighbor to two, three, or four existing clusters, it is not considered 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 1and 2-coalescences result in merging only two clusters; accordingly, most of the observed coalescences (about 95% according to our simulations on L = 2000 lattice) involve only two clusters while coalescence of three or four clusters is a rare event. The consecutive process of tree formation for a simple four-particle cluster is illustrated in Fig. 3. To construct the tree one needs to consider all consecutive coalescences that have formed the cluster, not

ð1Þ

2.2. Tree representation of clusters We represent each cluster 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 merge 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

Fig. 2. 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 simulations with L = 2000 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.

96

I. Zaliapin et al. / Tectonophysics 413 (2006) 93–107

Fig. 3. Tree representation of clusters: scheme. The dynamics is from left to right. At first step particle A is dropped onto the lattice and a oneparticle 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 corresponds to a five-node tree.

only its final shape. Therefore, it is clear that the same tree may correspond to clusters of different shape: Fig. 4a shows two 11-particle clusters that both correspond to the same tree shown in panel b. Accordingly, working with trees, we unavoidably narrow the information about the cluster population. Notice however that trees capture an excessively larger amount of information than mere cluster masses. Summing up, the time evolution of a cluster is necessary and sufficient to uniquely determine the corresponding tree, while the inverse is not true. Generally, clusters of different shape may correspond to the same tree (see Fig. 5); as well as clusters of the same final shape but with different formation history may correspond to different trees (see Fig. 4). The non-uniqueness of the tree representation of clusters does not affect our results, and it 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.

2.3. Horton–Strahler ranking 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 solution to this problem is given by the Horton–Strahler topological classification of ramified patterns (Horton, 1945; Strahler, 1957; Badii and Politi, 1997; Newman et al., 1997) illustrated in Fig. 4b: 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, ..., n z 2, merge together, a new cluster is formed with the rank (Badii and Politi, 1997):
if ri ¼ r1 8i ¼ 1; . . . ; n r1 þ 1; r¼ maxðri Þ; otherwise: The rank of a cluster is that of the root of the corresponding tree. It is also possible to consider an

Fig. 4. (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.

I. Zaliapin et al. / Tectonophysics 413 (2006) 93–107

97

Fig. 6. Mass–rank distribution observed on a 2000
2000 lattice at percolation. Dots—individual clusters, balls—average mass  r1 m r within a given rank. Line shows the relation mr ¼ 100:625 c4:2r1 :

Fig. 5. Non-uniqueness of tree representation. Three clusters with the same final shape but different history of formation are represented by different hierarchical trees. Particles have been dropped according to their alphabet marks; so first was the particle A, then B, etc. The Horton–Strahler ranks are shown next to the tree nodes.

alternative definition of ranks (Newman et al., 1997): 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 definitions coincide when only two clusters coalesce. The results reported in this paper are independent of a particular definition, since coalescence of more than two clusters (especially of high ranks) is a rare event.

proximately lognormal (not shown) with the mean given by Eq. (2) and a rank-independent standard deviation. The steady-state simulations of (Morein et al., 2004) suggest c = 4.325, which is reasonably close to c = 4.2 that we observe at percolation. The exponential relation of Eq. (2) happens to be valid over the entire time interval 0 b q b q c; the corresponding dynamics of c (q) is shown in Fig. 7. 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

3. Mass–rank relationship In this section we analyze the average mass m r of clusters with a given rank r; this is important for connecting various mass and rank scaling laws. The observed mass–rank distribution of clusters at percolation is shown in Fig. 6; it obeys the exponential relation mr ¼ 10cðr1Þ ¼ cr1 ;

ð2Þ

with c c 0.625, c = 10c c 4.2. Our simulations suggest that the mass distribution within a given rank is ap-

Fig. 7. Parameter c of the mass–rank relation mr = c r  1 as a function of time. At percolation c (q c) 6 4.2; the Euclidean limit of Gabrielov et al. (1999) corresponds to c = 3.25, it is reached at q c  q 6 0.14.

98

I. Zaliapin et al. / Tectonophysics 413 (2006) 93–107

not monotonous; it is accompanied by pronounced logperiodic 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 Newman et al., 1995), biased diffusion on random lattices (Stauffer and Sornette, 1998), diffusion-limited aggregation (DLA) (Sornette et al., 1996), and others. Log-periodic oscillations can be naturally explained by the Discrete Scale Invariance (DSI) (Sornette, 2004), 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 possesses a discrete scale symmetry: In order to superimpose its scaled image onto the original, one has to stretch it by the discrete factors 3n , n = 1, 2, ..., 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. Recall that the models of (Gabrielov et al., 1999; Morein et al., 2004) use the bfractal correctionQ e to the cluster shape; this correction affects the rate r ij of cluster coalescence: rij cejjij Li Lj ; where L i is the total boundary size of the clusters of rank i. The correction e can be expressed as (Gabrielov et al., 1999) 1 c1 ; e ¼ pffiffiffi c c2

a proper description of the time-dependent rank distribution, suggesting an original non-trivial character of a hierarchical description. In this section we summarize time-dependent laws for the mass- and rank-distributions (Stauffer and Aharony, 1994; Zaliapin et al., 2004), starting with a look at empirical observations. 4.1. Observed dynamics of cluster population The dynamics of the total number (n r d L 2) of the clusters of a given rank r is illustrated in Fig. 8 for r = 5, 6, 7. The population follows a characteristic bell-shaped trajectory with the peak prior to percolation. One does not observe a 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. Fig. 9 shows the population dynamics for the ranks 1 V r V 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 (Turcotte et al., 1999) that shows the dynamics of clusters with logarithmically binned masses. Figs. 10 and 11 illustrate the mass and rank distribution for a system with L = 2000 at q = 0.29 (dashed line), and at percolation (dash-dotted line). To smooth out statistical fluctuations, Fig. 10 shows the number of clusters with mass equal to or larger than m: P mV z m n m (q c). The cluster size distribution is prominently transformed during the model development: its downward bend at the early times is substituted at per-

which, together with results of Fig. 7, shows that in the percolation model the effective value of e decreases in time passing the Euclidean limit e = 1 (Gabrielov et al., 1999) at (q c  q) c 0.14 and approaching the steadystate bfractalQ e c 0.68 (Morein et al., 2004) at q = q c. The interval 2 b c V 4.2 observed during 0 b q V q c corresponds to l N e z 0.68. 4. Cluster size distributions Time-dependent rank distribution of clusters provides essential information about the model dynamics. Also it presents a direct counterpart of the Gutenberg– Richter law for magnitude distribution in observed seismicity. The study (Zaliapin et al., 2004) demonstrated that simple substitution of the mass–rank relation (2) into well-known laws for the dynamics of cluster masses (Stauffer and Aharony, 1994) does not provide

Fig. 8. Dynamics of population n r d L 2 of a given rank, r = 6, 7, 8 for L = 2000. Moment of percolation is depicted by a vertical dashed line.

I. Zaliapin et al. / Tectonophysics 413 (2006) 93–107

Fig. 9. Dynamics of population n r of a given rank, 1 V r V 11 in semilogarithmic scale. Moment of percolation is depicted by a vertical dashed line. (Cf. Fig. 7 from Turcotte et al., 1999).

colation with an almost linear (in bi- or semi-log scale) shape with a slight upward bend for the largest clusters. 4.2. Time-dependent laws for cluster size distribution It is known that the time-dependent mass distribution is described by Stauffer’s two-exponent scaling law (Stauffer, 1975; Stauffer and Aharony, 1994; Margolina et al., 1984): s

nm ðqÞfm f0 ð zÞ;

z ¼ ðqc  qÞm

1=2

þ z0 ;

ð3Þ

where n m (q) is the number of clusters of mass m per lattice site at instant q. The Fisher exponent s = 187 / 91 c 2.05 is universal for 2D systems (Fisher, 1967;

99

Fig. 11. Rank distribution of clusters observed for 2000
2000 lattice at percolation q = q c (dash-dotted line), q = 0.29 (dashed line), and averaged over the percolation cycle 0 b q b q c (solid line). For comparison, all curves are normalized to unity at r = 1.

Stauffer and Aharony, 1994); the function f 0 has a bellshaped form with maximum to the left of percolation; and the shift z 0 is independent of m. Considered as a function of m, the two-exponent scaling explains the power law mass distribution at percolation (with q 0 = f 0 (z 0)): nm ðqc Þfq0 ms ;

ð4Þ

it also explains the downward bend for q b q c, clearly observed in Fig. 10 (dashed line). At the same time, as a function of q it describes the bell-shaped dynamics of clusters with given mass m. Eq. (4) suggests the slope s  1 c1.05 for the cumulative mass distribution, while the observed slope 0.96 is somewhat less than that (see Fig. 10). This is due to the impact of two concurrent phenomena: so-called bdeviation from scalingQ at small m (Hoshen et al., 1979) and finite-size effects at large m (Margolina et al., 1984; Hoshen et al., 1979) (see Section 4). When we switch to the rank description, the twoexponential scaling is no longer sufficient. As shown in (Zaliapin et al., 2004), the number n r (q) of clusters of rank r at instant q is obeying the three-exponent scaling nr ðqÞfg0 ð zÞ10br ; hðrÞ ¼ a1 10r1 r ;

z ¼ ðqc  qÞhðrÞ þ zV0 ; zV0 ðrÞ ¼ a2 10r2 r

ð6Þ

with Fig. 10. Mass distribution of clusters observed on a 2000
2000 lattice at percolation q = q c (dash-dotted line), q = 0.48 (dashed line), and averaged over 0 b q b q c (solid line). To smooth Pout statistical fluctuations we show the cumulative distribution: ~ mVNm n m d . For comparison, all curves are normalized to unity at m = 1.

r1 c0:23;

r2 c0:03;

a1 c1:54;

ð5Þ

a2 c1:43

and function g 0 being approximately Gaussian:   g0 ð zÞ ¼ exp  z2 =2 :

100

I. Zaliapin et al. / Tectonophysics 413 (2006) 93–107

The three-exponent scaling (5) is consistent with an exponential approximation to the rank distribution at percolation (with P 0 = g 0 (zV0)) nr ðqc Þfp0 10br ;

bc0:62:

ð7Þ

Possible deviations from the pure exponential law at q b q c (clearly observed in Fig. 11) and dynamics of a given rank (see Figs. 8 and 9) are described by specific form of the functions g 0 (d ), h (d ), and zV0 (d ). The function g 0 (z) is shown in Fig. 12 where different symbols depict clusters of different ranks. The collapse is obvious, confirming the validity of the three-exponent scaling (5) and (6). In the Stauffer’s scaling (3) for cluster masses, the time renormalization (q c  q)m 1 / 2 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 (5) and (6) for ranks the time renormalization (q c  q)10r 1r collapses the dynamics of rank r clusters onto the master curve g 0 (z  zV0), although the shift now is rank dependent and is given by 10r2r. 4.3. Averaged scaling 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. (Daley and Vere-Jones, 2003). In such situations one uses the averaged measurements.

Fig. 12. Three-exponent scaling for rank dynamics. The master Gaussians g 0 (z) for different ranks collapse when using the renormalization given by Eqs. (5) and (6).

For instance, the famed Gutenberg–Richter law (Gutenberg and Richter, 1954; Turcotte, 1997; Ben-Zion, 2003) that gives exponential approximation to the size distribution of earthquakes (via their magnitudes) is valid only after appropriate averaging over a suitable spatio-temporal domain. This explains the importance of the question: How do the distributions of Eq. (3) and (5) change after temporal averaging? The integration of (3) over the interval 0 V q V q c gives (Zaliapin et al., 2004): ^ nm~m5=2 :

ð8Þ

The validity of (8) is confirmed by the observed averaged mass distribution shown by the solid line in Fig. 10: it retains the power-law form (4) of the distribution at percolation while the slope is increased by 1 / 2 due to averaging. Similarly, for ranks one obtains (Zaliapin et al., 2004): ^ nr ~10ra ;

ac0:87:

ð9Þ

Again, the averaged rank distribution retains the exponential form (7) of the distribution at percolation; while its index has increased due to averaging (see Fig. 11). 4.4. Correction to simple scaling Due to finiteness of the lattice, the results of previous sections require some corrections to exactly match the simulated rank distributions. The pure power and exponential laws in Figs. 10 and 11 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 bdeviations from scalingQ (Hoshen et al., 1979): it can be shown analytically that the small clusters do not yet obey the general scaling law of Eqs. (4) and (7) which holds only for large enough masses (ranks). The upward bend at large clusters is due to finite-size effects (Hoshen et al., 1979; Margolina et al., 1984): each large cluster that reaches outside the lattice boundary is bseenQ 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. (1979) and Margolina et al. (1984).

I. Zaliapin et al. / Tectonophysics 413 (2006) 93–107

In case of the mass distribution, the corrections to scaling are given by (Margolina et al., 1984):  nm ðqc Þgms q0 þ q1 mX þ qL m1=D L1 ; ð10Þ where X c 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 the rank distribution, the bdeviations from scalingQ at lower clusters are only observed for r = 1; while the finite-size effects at large clusters are clearly present for many ranks. Zaliapin et al. (2004) proposed the following correction to scaling for the rank distribution:   nr ðqc Þg10br p0 þ pL 10dr L1 ;rN1: ð11Þ with d¼

1 log cc0:33 D 10

and showed that the observed ranks 4 V r V 9 nicely follow the predicted scaling. Importantly, the corrections to scaling (11) act at all cluster sizes, so they cannot 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. (2004) (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 Figs. 10 and 11 and is also seen in the analysis of Turcotte et al. (1999) (their Fig. 9).

101

ditionally and effectively used in earthquake prediction research. We emphasize that studying the percolation cluster may differ from studying the largest cluster on the grid. Indeed, for large lattices the percolation cluster is almost surely is the largest one at the moment of percolation; this equivalence though may not hold at the earlier stages of the model dynamics as well as for intermediate size lattices. Below we use the term bpercolation clusterQ to denote the cluster that will percolate first in a given model run. 5.1. Mass dynamics of percolation cluster The time-dependent mass m p (q) of the percolation cluster, averaged over 1000 runs of the model on L = 2000 grid, is shown in Fig. 13. One can easily distinguish the following three major stages of the mass development: Intermediate dynamics, 0.3 N (q c  q) N 0.01, is well described by the scaling law mp c5:0
ðqc  qÞ2:05 :

ð12Þ

However, at early times, (q c  q) N 0.3, as well as in the immediate vicinity of percolation, (q c  q) N 0.01, the mass dynamics deviates downward from the scaling of Eq. (12). At early times this is due to the bdeviations from scalingQ (Hoshen et al., 1979), while at later times to the finite size effects (Hoshen et al., 1979; Margolina et al., 1984) (see Section 4). The downward deviation from the power scaling at later times seems a nice signal of the approaching percolation. However such a signal might be undetectable in a particular model run, as illustrated in the insert of Fig. 13: While the overall scaling with exponent 2.05 is clearly observed,

5. Properties of percolation cluster This section analyzes the growth of the percolation cluster; that is, we consider the time-dependent properties of the cluster that will first span the lattice (unavoidably using the information from the bfutureQ). Such properties, being compared to the average statistics over the entire cluster population, might be useful for general description of the model as well as for improving the observer’s knowledge about the moment of the approaching percolation. Studying the percolation cluster separately from the rest of the model can be considered an analog to concentrating on the region surrounding the epicenter of a large regional earthquake versus the entire region; the latter approach being tra-

Fig. 13. Dynamics of the mass of the cluster which will percolate first. Main panel shows the average over 1000 runs of a model with L = 2000, insert shows the dynamics of a single model run.

102

I. Zaliapin et al. / Tectonophysics 413 (2006) 93–107

the downward deviation prior to percolation is not seen for this particular realization. This problem is very similar to the one typically faced in earthquake prediction studies: many premonitory seismicity patterns (the Benioff strain release (Jaume and Sykes, 1999) being the most popular one) are easily detected by the stacking data from many large earthquakes, but almost unseen in the observed seismicity prior to a single earthquake. It should be noted though that the direct analog of a power law Benioff strain release—a power law mass increase of the percolation cluster—seems to be nicely detectable even in a single model run (see insert in Fig. 13), and may serve as a basis for predicting the onset of percolation. We leave a detailed treatment of such a prediction technique for a further study, and proceed now with the dynamics of the rank of the percolation cluster. 5.2. Rank dynamics of percolation cluster We argued above that the ranks are more feasible for observations than masses. This advocates the study of rank dynamics for the percolation cluster. Fig. 14 shows the time-dependent rank r p(q) of the percolation cluster, averaged over 1000 runs of the model on an L = 2000 grid. The following exponential increase is clearly observed at intermediate times to percolation, 0.3 N (q c  q) N 0.01: rp c2:2  3:35
log10 ðqc  qÞ:

ð13Þ

As in the mass dynamics, the downward deviations from this law at early times, (q c  q) N 0.3, and close to percolation, (q c  q) N 0.01 can be explained by bdevia-

Fig. 14. Dynamics of the rank of the cluster which will percolate first; average over 1000 runs of a model with L = 2000.

tions from scalingQ and finite size effects (see Section 1). Again, the delay in creation of the consecutive rank may serve as a signal of approaching percolation, but this delay would be hard (if possible at all) to detect in a single model run. Striking is the following observation: The building of the percolation rank is monotonous, it lacks the characteristic log-periodic oscillations observed for the process of rank creation on the entire grid (see Fig. 7 and its discussion). 5.3. Mass–rank relation for percolation cluster Having established the scaling laws for mass (Eq. (12)) and rank (Eq. (13)) of the percolation cluster, we now proceed with connecting these laws. One obviously may rewrite them as  m 1=2:05 p ðqc  qÞc ; ðqc  qÞc10ðrp 2:2Þ=3:35 ; 5:0 which leads to  m 1=2:05 p c10ðrp 2:2Þ=3:35 5:0 or

 2:05 2:05
2:2  log10 5:0  log10 mp crp 3:35 3:35 c0:61
rp  0:65: The index 0.61 is fairly close to that of the mass– rank relation (2) (with c c 0.625) estimated for the entire cluster population at percolation (see Fig. 6). This prompts one to conjecture that the index of the mass–rank scaling for the percolation cluster matches that of the total cluster population at percolation. This observation, schematically illustrated in Fig. 15, depicts a very special character of the percolation cluster: at any time instant it possesses the information about the index of the mass–rank scaling law that will be observed for the entire cluster population at percolation. Thus, comparison of the mass–rank scalings for the entire cluster population with that of the percolation cluster may serve as another precursor for the percolation onset. Indeed, we do not know in advance which cluster will percolate first, so one may use instead the largest observed cluster. Detailed analysis of the power and usefulness of such a premonitory pattern is beyond the scope of the present study. We also observe from Fig. 15 that the mass m p of the cluster that will percolate first is always less than the average mass of clusters of the same rank at percolation. This observation is supported by Fig. 16 that shows the dynamics of the average mass m r of clusters

I. Zaliapin et al. / Tectonophysics 413 (2006) 93–107

Fig. 15. Mass–rank relationship for the entire cluster population at percolation (top solid line), entire cluster population prior to percolation (dashed line), and the cluster that will percolate first (bottom solid line).

103

Fig. 17. Rank distribution of subclusters within the cluster that will percolate first. For comparison, all curves are normalized to unity at r = 1 and arbitrarily shifted downward.

of a given rank. Strikingly, the average mass continues to increase when both the total number n r of clusters and their total mass M r significantly decreased having passed their respective peak values. In other words, the bfatQ members of a rank population still survive when their bthinQ companions already have merged together turning into clusters of a higher rank.

steps. Thus, it is natural to consider the rank distribution of subclusters within the percolation cluster at different stages of its development. This distribution for four time instants: q = 0.29, 0.5, 0.54, and at percolation is shown in Fig. 17. The rank distribution within the percolation cluster is always approximately exponential (modulo the largest ranks that deviate significantly from the overall pattern):

5.4. Rank distribution within percolation cluster

f m r ~100:3r ;

Recall that each cluster in our model is an aggregate of smaller ones that merged together at previous time

ð14Þ

˜ r, is the number of clusters of rank r within the where m percolation cluster. Notable is the upward bend that reflects the bcriticalQ character of the percolation cluster (see Section 7 for details). Importantly, the distribution is stable in time; there is no transformation from convex to almost linear shape observed for the entire cluster population (see Section 4, Fig. 11). The index of the exponential law (14) is twice less than that of the rank distribution for the entire cluster population (see Section 4, Eq. (7), Fig. 11). Thus the relative number of larger clusters within the percolation one is higher comparing to that for the entire cluster population. 6. Tokunaga scaling for side branching

Fig. 16. 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.

In this section we show how the index c of the mass–rank relationship (2) is expressed via the parameters of another important scaling law (Tokunaga scaling (Tokunaga, 1978)) that describes self-similar side branching of clusters. First, following (Newman et al., 1997; Tokunaga, 1978) we define the branching ratio Tij for a given cluster (tree) as the number N ij of its subclusters (nodes) of rank i that joined subcluster

104

I. Zaliapin et al. / Tectonophysics 413 (2006) 93–107

(node) of rank j, averaged over subclusters (nodes) of rank j: Tij ¼

Nij : Nj

Many natural hierarchies are shown to exhibit selfsimilar side branching, obeying the Tokunaga scaling constraint (Newman et al., 1997; Tokunaga, 1978): Tiiþk ¼ Tk ¼ s0 s

k1

ð15Þ

:

Particularly, it was predicted by Gabrielov et al. (1999) and later confirmed by simulations (Morein et al., 2004) that clusters produced by aggregation dynamics obey the Tokunaga scaling (15) asymptotically in k. Moreover, Gabrielov et al. (1999) established a relationship between parameters s 0,s of the Tokunaga scaling (15) and the mass–rank distribution (2). Specifically, if one assumes clusters of a regular non-fractal shape (Euclidean limit of the model), then s ¼ 1=s0 c1:80193774;

s0 c0:55495813; c¼

1=s20 c3:24697602:

ð16Þ

Observed systems though are expected to produce fractal clusters (non-Euclidean case), which might change the values of the parameters and their interplay; for example, the steady-state simulations of Morein et al. (2004) suggest sc3:0253;

s0 c0:6993;

cc4:325

ð17Þ

The general connection between s 0, s and c can be established by noticing that the Tokunaga scaling uniquely determines the mass of rank r clusters. To express the average mass m r via the parameters s 0, s we notice 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: m1 ¼ 1 m2 ¼ ð2m1 þ PÞ þ T12 ðm1 þ PÞ m3 ¼ ð2m2 þ 1Þ þ T23 ðm2 þ 1Þ þ T13 ðm1 þ PÞ ... k1 X mk ¼ð2mk1þ1Þþ Tkik ðmki þ1Þ  ð1 PÞT1k ;kz3: i¼1

ð18Þ 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 N 2 can only join other clusters using a one-particle connector. Assuming the Tokunaga scaling (15) we rewrite Eq. (18) for k z 3 as k 1 X Ti ðmki þ 1Þ  ð1  PÞTk1 ; mk ¼ ð2mk1 þ 1Þ þ i¼1

and using the mass–rank relation (2) with c N 1 obtain ck1 ¼ 2ck2 þ 1 " ¼ ck2 2 þ


k1 X

  s0 si1 cki1 þ 1 ð1  PÞs0 sk1

i¼1

1 ck2

þ s0

1  ðs=cÞk1 s0 þ k2 1  s=c c #

sk1  1  ð1  PÞs0 ðs=cÞk2 s1

It is easily checked that this equation has a solution only if c N s; thus s/c b 1 and for large k then follows   s0 ck1 ¼ ck2 2 þ 1  s=c leading to the final equation c2  cð2 þ s þ s0 Þ þ 2s ¼ 0: A simple analysis shows that this equation always has at least one solution for s 0 N 0, and the only solution that ensures c N s is given by qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 þ s þ s0 þ ð2 þ s þ s0 Þ2  8s c¼ ; ð19Þ 2 while c N 1 implies s N 1 + s 0. For the Euclidean values of s 0, s (16) the above formula gives c (s 0, s) = 3.24697960, which is remarkably close (6 digits) to the result of (Gabrielov et al., 1999). Furthermore, the non-Euclidean values of (17) solve Eq. (19) exactly. We found it quite convincing that our complimentary set of assumptions used to derive (19) lead to the same numerical results as an analytical study (Gabrielov et al., 1999) and simulations of (Morein et al., 2004). This suggests an internal stability of a system that follows an aggregation dynamics with respect to the details of its particular realization. The relationship (19) provides a striking constraint for the scaling of aggregation processes, and may be useful in working with observed systems. 7. Discussion This work continues a line of research initiated recently by Turcotte et al. (1999), Gabrielov et al.

I. Zaliapin et al. / Tectonophysics 413 (2006) 93–107

(1999), Morein et al. (2004), and Zaliapin et al. (2004), who have introduced and studied hierarchical inverse cascade description of aggregation (coagulation) phenomena. Notably, the inverse cascade approach is shown to explain self-organized critical behavior (Turcotte et al., 1999), while the aggregation phenomena are proven important for evolution of many natural hazardous processes such as earthquakes, landslides, and forest fires (Turcotte et al., 2002; Morein et al., 2004); they are also commonly employed for description and modeling of material fracture. Following (Gabrielov et al., 1999) we described clusters in 2D site-percolation model on a square lattice 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 (Gabrielov et al., 1999; Morein et al., 2004) that have developed mean-field steady-state approximations to the system. The most important results of percolation theory are traditionally given in terms of cluster masses (Stauffer and Aharony, 1994). Accordingly, in order to bridge the classical results to the novel rank description it is vitally important to establish the mass– rank scaling laws, which was done in (Gabrielov et al., 1999; Morein et al., 2004). In this study we have analytically expressed the index of the mass–rank exponential relationship (2) via the parameters of the Tokunaga side-branching constraint (15). The predictions and simulations of Gabrielov et al. (1999) and Morein et al. (2004) are found to be in a perfect agreement with this relationship, suggesting the universality of the essential characteristics of the aggregation process, independently of its particular (physical or analytical) realization. We studied in detail the transition of the system from earlier stages to the vicinity of percolation and reported several characteristic phenomena observed as q Y q c. They include transformation of the cluster size distribution not unlike that observed in seismicity, steel samples fracture, and previous models of hierarchical fractures (Narkunskaya and Shnirman, 1990; Rotwain et al., 1997; Gabrielov et al., 2000a; Zaliapin et al., 2003). 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

105

described by intermediate depth hierarchies (in other words they have an 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. (2004) show that finite-size effects are not negligible even for large L. We have studied the growth of the percolation cluster tracing its evolution from the earliest times to the onset of percolation and established scaling laws for mass and rank of the percolation cluster. We also studied the frequency-size distribution of subclusters within the percolation cluster. We have reported several features that may be used for prediction of the onset of percolation: power-law increase of the mass of the percolation cluster (this property is known from the classical studies (Stauffer and Aharony, 1994)); exponential increase of the rank of the percolation cluster; downward deviations of the mass- and rankscalings prior to percolation; and convergence of the index of the mass–rank distribution for the entire cluster population to that within the percolation cluster. Indeed, if we may count precisely the occupied sites, then the onset of percolation is determined by the critical time q c, and there is no need in additional premonitory patterns. At the same time, it is not uncommon in prediction research that premonitory patterns whose emergence is obvious in a simple model can be effectively applied for prediction of critical events in more complex systems as well as in observations (e.g., see Gabrielov et al., 2000b; Narkunskaya and Shnirman, 1990; Narkunskaya and Shnirman, 1994; Rotwain et al., 1997; Shebalin et al., 2000; Zaliapin et al., 2003), so the further analysis of the above patterns seems promising. Our closing remark is on the index s of cluster mass distribution at percolation (Eq. (4)). The studies of Gabrielov et al. (1999) and Morein et al. (2004) predict s = 2; which slightly deviates from the well established theoretical value of the Fisher exponent s = 187 / 91 c 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 (Gabrielov et al., 1999; Morein et al., 2004). In fact, this implication is not true. To validate the approach of (Gabrielov et al., 1999; Morein et al., 2004) 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 (Ziff et al., 1999; Cardy and Ziff,

106

I. Zaliapin et al. / Tectonophysics 413 (2006) 93–107

2003) have demonstrated that when we bcharacterize the size distribution of clusters in a way that circumvents the site-level descriptionQ considering any bmacroscopic measure of the length scale of the clusterQ, the exponent of the corresponding scaling law becomes 2, universally for all 2D systems. An example of a bmacroscopic measureQ is the linear size in an arbitrary direction, the radius of gyration, the diameter of the covering disk, etc. Clearly, the modeling of (Gabrielov et al., 1999; Morein et al., 2004) deals with such a macroscopic measure of cluster size, and hence predicts the correct scaling exponent. Acknowledgments We are sincerely grateful to Bill Newman for numerous focused discussions that have helped significantly in organizing and presenting the results. We thank Kristy Tiampo and two anonymous reviewers for their comments that helped to improve the paper; and Vladimir Keilis-Borok, Gleb Morein, and Donald Turcotte for their continuous interest in the work. This study was partly supported by NSF, grant ATM 0327558. References Aki, K., 1995. Earthquake prediction: societal implications. Rev. Geophys., 243 – 247. Allegre, C.J., LeMouel, J.L., Provost, A., 1982. Scaling rules in rock fracture and possible implications for earthquake prediction. Nature 5861, 47 – 49. Badii, R., Politi, A., 1997. Complexity: Hierarchical Structures and Scaling in Physics. Cambridge University Press. 318 pp. Bak, P., Tang, C., Wiesenfeld, K., 1998. Self-organized criticality. Phys. Rev. A, 364 – 374. Ben-Zion, Y., 2003. Appendix 2, Key Formulas in Earthquake Seismology. International Handbook of Earthquake and Engineering Seismology: Part B. Academic Press, pp. 1857 – 1875. 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. 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. 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. 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. Stud. Appl. Math. 3, 163 – 204. Daley, D.J., Vere-Jones, D., 2003. An Introduction to the Theory of Point Processes, 2nd ed., vol. I,. Springer, New-York, p. 469. Fisher, M.E., 1967. Physics, 255. Gabrielov, A., Newman, W.I., Turcotte, D.L., 1999. An exactly soluble hierarchical clustering model: inverse cascades, self-similarity, and scaling. Phys. Rev. E, 5293 – 5300.

Gabrielov, A., Keilis-Borok, V., Zaliapin, I., Newman, W.I., 2000a. Critical transitions in colliding cascades. Phys. Rev. E, 237 – 249. Gabrielov, A.M., Zaliapin, I.V., Keilis-Borok, V.I., Newman, W.I., 2000b. Colliding cascades as a model for earthquake prediction. Geophys. J. Int., 427 – 437. Gutenberg, B., Richter, C.F., 1954. Seismicity of the Earth and Associated Phenomena, 2nd ed. Princeton University Press, Princeton. Horton, R.E., 1945. Erosional development of streams and their drainage basins: hydrophysical approach to quantitative morphology. Geol. Soc. Am. Bull., 275 – 370. Hoshen, J., Stauffer, D., Bishop, G.H., Harrison, R.J., Quinn, G.D., 1979. Monte Carlo experiments on cluster size distribution in percolation. J. Phys. A, 1285 – 1303. Jaume, S.C., Sykes, L.R., 1999. Evolving towards a critical point: a review of accelerating seismic moment/energy release prior to large and great earthquakes. Pure Appl. Geophys. 2–4, 279 – 305. Kadanoff, L.P., 2000. Statistical Physics: Statics, Dynamics, and Renormalization. World Scientific Publishing, Singapore. Keilis-Borok, V.I., 1996. Intermediate-term earthquake prediction. Proc. Natl. Acad. Sci. U. S. A., 37 – 3755. Keilis-Borok, V., 2002. Earthquake prediction: state-of-the-art and emerging possibilities. Ann. Rev. Earth Planet. Sci., 1 – 33. Keilis-Borok, V.I., Soloviev, A.A. (Eds.), 2003, Nonlinear Dynamics of the Lithosphere and Earthquake Prediction. Springer, Heidelberg. Knopoff, L., Newman, W.I., 1983. Crack fusion as a model for repetitive seismicity. Pure Appl. Geophys. 3, 495 – 510. Leyvraz, F., 2003. Scaling theory and exactly solved models in the kinetics of irreversible aggregation. Phys. Rep. 2–3, 95 – 212. Margolina, A., Nakanishi, H., Stauffer, D., Stanley, H.E., 1984. Monte Carlo and series study of corrections to scaling in twodimensional percolation. J. Phys. A, 1683 – 1701. Morein, G., Newman, W.I., Turcotte, D.L., Gabrielov, A., 2004. An Inverse Cascade Model for Self-Organized Complexity and Natural Hazards, p. 9. Preprint. Narkunskaya, G.S., Shnirman, M.G., 1990. Hierarchical model of defect development and seismicity. Phys. Earth. Planet. Inter., 29 – 35. Narkunskaya, G.S., Shnirman, M.G., 1994. An algorithm of earthquake prediction. Computational Seismology and Geodynamics. AGU, Washington, D.C, pp. 20 – 24. Newman, W.I., Gabrielov, A., 1991. Failure of hierarchical distributions of fiber bundles. I. Int. J. Fract., 1 – 14. Newman, W.I., Knopoff, L., 1982. Crack fusion dynamics—a model for large earthquakes. Geophys. Res. Lett. 7, 735 – 738. Newman, W.I., Knopoff, L., 1983. A model for repetitive cycles of large earthquakes. Geophys. Res. Lett. 4, 305 – 308. Newman, W.I., Knopoff, L., 1990. Scale-invariance in brittle-fracture and the dynamics of crack fusion. Int. J. Fract. 1, 19 – 24. Newman, M.E.J., Ziff, R.M., 2001. Fast Monte Carlo algorithm for site or bond percolation. Phys. Rev. E, 016706. Newman, W.I., Turcotte, D.L., Gabrielov, A., 1995. Log-periodic behavior of a hierarchical failure model with applications to precursory seismic activation. Phys. Rev. E, 4827 – 4835. Newman, W.I., Turcotte, D.L., Gabrielov, A., 1997. Fractal trees with side branching. Fractals, 603 – 614. Rotwain, I., Keilis-Borok, V., Botvina, L., 1997. Premonitory transformation of steel fracturing and seismicity. Phys. Earth Planet. Inter., 61 – 71. Rundle, J., Turcotte, D., Klein, W. (Eds.), 2002. Geocomplexity and the Physics of Earthquakes. American Geophysical Union, Washington, DC.

I. Zaliapin et al. / Tectonophysics 413 (2006) 93–107 Shebalin, P., Zaliapin, I., Keilis-Borok, V., 2000. Premonitory raise of the earthquakes correlation range: lesser Antilles. Phys. Earth Planet. Int. 3–4, 241 – 249. Shnirman, M.G., Blanter, E.M., 2001. Criticality in a Dynamic Mixed System. Phys. Rev. E 5 Art. No. 056123: Part 2. Sornette, D., 2004. Critical Phenomena in Natural Sciences, 2nd ed. Springer-Verlag, Heidelberg, p. 528. Sornette, D., Johansen, A., Arneodo, A., Muzy, J.-F., Sauleur, H., 1996. Complex fractal dimensions describe the internal hierarchical structure of DLA. Phys. Rev.Lett., 251 – 254. Stauffer, D., 1975. Violation of dynamical scaling for randomly dilute Ising ferromagnets near percolation threshold. Phys. Rev. Lett. 6, 394 – 397. Stauffer, D., Aharony, A., 1994. Introduction to Percolation Theory, 2nd ed. Taylor & Francis. Stauffer, D., Sornette, D., 1998. Log-periodic oscillations for biased diffusion on random lattice. Physica, A 3–4, 271 – 277. Strahler, A.N., 1957. Quantitative analysis of watershed geomorphology. Trans. Am. Geophys. Un., 913 – 920. Tokunaga, E., 1978. Consideration on the composition of drainage networks and their evolution. Geogr. Rep. Tokyo Metrop. Univ., 1 – 27. Toroczkai, Z., 2001. Topological classification of the Horton–Strahler index on binary trees. Phys. Rev. E, 016130.

107

Turcotte, D.L., 1991. Earthquake Prediction. Ann. Rev. Earth Planet. Sci., 263 – 281. Turcotte, D.L., 1997. Fractals and Chaos in Geology and Geophysics, 2nd edition. Cambridge University Press. 398 pp. Turcotte, D.L., 2001. Self-organized criticality: does it have anything to do with criticality and is it useful? Nonlinear Proc. Geophys., 193 – 196. Turcotte, D.L., Pelletier, J.D., Newman, W.I., 1998. Networks with side branching in biology. J. Theor. Biol. 4, 577 – 592. Turcotte, D.L., Malamud, B.D., Morein, G., Newman, W.I., 1999. An inverse cascade model for self-organized critical behavior. Physica, A, 629 – 643. Turcotte, D.L., Malamud, B.D., Guzzetti, F., Reichenbach, P., 2002. Self-organization, the cascade model, and natural hazards. Proc. Natl. Acad. Sci., 2530 – 2537. 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. 3–4, 839 – 861. Zaliapin I., Wong, H., Gabrielov, A., 2004. Inverse cascade in percolation model: hierarchical description of time-dependent scaling. To appear in Phys. Rev. E. Ziff, R.M., Lorenz, C.D., Kleban, P., 1999. Shape-dependent universality in percolation. Physica, A 1–4, 17 – 26.