Creator

Date

2000-08-29

Description

Comment: 12 pages, 13 figures, revtex+epsfig style, to appear in Phys. Rev. E

(Nov. 2000)

(Nov. 2000)

The existence and stability of the universality class associated to local

minimal energy landscapes is investigated. Using extensive numerical

simulations, we first study the dependence on a parameter $\gamma$ of a partial

differential equation which was proposed to describe the evolution of a rugged

landscape toward a local minimum of the dissipated energy. We then compare the

results with those obtained by an evolution scheme based on a variational

principle (the optimal channel networks). It is found that both models yield

qualitatively similar river patterns and similar dependence on $\gamma$. The

aggregation mechanism is however strongly dependent on the value of $\gamma$. A

careful analysis suggests that scaling behaviors may weakly depend both on

$\gamma$ and on initial condition, but in all cases it is within observational

data predictions. Consequences of our results

minimal energy landscapes is investigated. Using extensive numerical

simulations, we first study the dependence on a parameter $\gamma$ of a partial

differential equation which was proposed to describe the evolution of a rugged

landscape toward a local minimum of the dissipated energy. We then compare the

results with those obtained by an evolution scheme based on a variational

principle (the optimal channel networks). It is found that both models yield

qualitatively similar river patterns and similar dependence on $\gamma$. The

aggregation mechanism is however strongly dependent on the value of $\gamma$. A

careful analysis suggests that scaling behaviors may weakly depend both on

$\gamma$ and on initial condition, but in all cases it is within observational

data predictions. Consequences of our results

Type

Database

Link to record

Show preview

Hide preview

ar
X

iv :c

on d-

m at

/0 00

84 24

v1 [

co nd

-m at.

sta t-m

ec h]

2 9 A

ug 20

00 Local minimal energy landscapes in river networks

Achille Giacometti INFM Unita´ di Venezia, Dipartimento di Scienze Ambientali, Universita´ di Venezia

Calle Larga Santa Marta DD2137, I-30123 Venezia-Italy

(February 1, 2008)

The existence and stability of the universality class as- sociated to local minimal energy landscapes is investigated. Using extensive numerical simulations, we first study the de- pendence on a parameter γ of a partial differential equation which was proposed to describe the evolution of a rugged land- scape toward a local minimum of the dissipated energy. We then compare the results with those obtained by an evolution scheme based on a variational principle (the optimal channel networks). It is found that both models yield qualitatively similar river patterns and similar dependence on γ. The ag- gregation mechanism is however strongly dependent on the value of γ. A careful analysis suggests that scaling behaviors may weakly depend both on γ and on initial condition, but in all cases it is within observational data predictions. Con- sequences of our results are finally discussed and the most plausible scenario is presented.

92.40.Fb,64.60.Ht,05.60.+w

I. INTRODUCTION

Understanding the development of a landscape in terms of fundamental mechanical principles is a formidable task [1]. In spite of this high complexity, re- cent theoretical studies have resulted in considerable pro- gresses by considering the issue from a viewpoint anal- ogous to the one taken in conventional critical phenom- ena [2–8] where simple models are exploited to identify universality classes. The main idea of this approach is indeed to focus on few fundamental ingredients which, in the spirit of critical phenomena, are expected to provide a reasonable description of the large-scales properties of the network. The remarkable properties of river networks have been

known for some time and they can be condensed in few phenomenological scaling laws which have been con- firmed to hold in the observational data [9]. While these laws do not explain the underlying physical mech- anisms, they nevertheless provide guidelines for their search. Hence any newly proposed model for river net- works ought to be tested against these laws. In the lan- guage of critical phenomena, those scaling laws can be used to derive critical exponents and thus discriminate among different universality classes. Among such laws a fundamental role in the physics of the erosion of a land- scape due to the flow of water over it, is played by the slope-area law (see Ref. [1] for a review).

A few Langevin-like equations have been recently proposed for describing the evolution of the landscape under the effect of erosional processes. In Ref. [5], reparametrization invariance arguments [10] were used to derive a dynamical equation which yields the slope- area law as a stationary state. The same equation was obtained by Somafai and Sander [7] using Landau ar- guments. Other proposals have also been advanced [6], [11–13]. The equation proposed in Ref. [5], was studied both

analytically in 1+1 dimensions (one spatial and one tem- poral) and numerically using a self-consistent (SC) solu- tion in 2+1 dimensions, and it was shown to predict a fairly reasonable stationary state quite different from the starting network acting as initial condition. This type of analysis hinges on the observation that the network rep- resenting the river has an evolution mirroring that of the evolution profile. Other methods have also been devised which directly address the topological properties of the network itself, the best known being the optimal channel networks (OCN) [14]. This is a lattice model where a functional describing the dissipated energy is minimized in order to find the optimal configuration and it is based on the idea that, presumably, the erosional process tak- ing place on a landscape is driven by a strive for opti- mality. A simulating annealing procedure [15] has been implemented to this aim both at zero [2] and finite tem- perature [16]. While the latter is aimed to find the ab- solute minimum, the former is expected to display local minima only. Using exact bounds and finite-size scaling, Maritan et al ( [17,18]) showed that the absolute mini- mum belongs to the mean-field universality class, that, in turn, means that the corresponding network has a highly symmetric pattern with small rivers draining into bigger rivers in a predictable way (this network is akin to the Peano basin [19]). However this minimum is not eas- ily reachable in the space of all configurations, and one is then led to suspect that real rivers are better described by configurations related to local minima. We further note that stationary states of the aforementioned dynamical equations, are also expected to be associated with local minima for the same reason. In view of the above discussion, it is apparent that a

more complete description of the stability and the scaling behavior associated with these local minima would be de- sirable. The present work is an attempt in this direction. Our aims are threefold. Our first goal is the full char- acterization of the possible universality class associated with local minima. Using the SC numerical procedure

1

envisaged in Ref. [5], we extend previous numerical re- sults both in size and statistics and we use a finite size procedure for a more accurate estimate of the exponents. We then compare these exponents with those calculated using a zero-temperature simulation in the OCN frame- work with similar sizes and statistics. Our second goal is the generalization of the dynamical equation to include a tunable parameter which turns out to be related to the 0 ≤ γ ≤ 1 parameter considered in OCNs and which is responsible of the aggregation mechanism. In Ref. [18] it was shown that the absolute minimum is insensitive to a variation of γ in the interval 1/2 ≤ γ < 1 [20]. We find that although the final patterns display marked differences as a function of γ, critical exponents show a much smaller discrepancy, which our results indicate to be marginal (at the edge of error bars), both in the SC and OCN case. Finally we test the stability of our results against a variation in the initial conditions. We find that even in this case a marginal difference appears in the effective critical exponents. The paper is organized as follows. In Sec.II, we re-

view the dynamical equation whereas in Sec.III scaling laws and critical exponents are briefly recalled. Sec.IV contains the definition of OCN and our results are pre- sented in Sec.V. Finally Sec.VI contains some concluding remarks and future perspectives.

II. THE DYNAMICAL EQUATION

A rather general dynamical equation consistent with general principles, and capturing the physics of erosional processes occurring in river basin, is [5], [7]:

∂h(x, t)

∂t = u+ ν∇2h(x, t)−Qα(x, t)

[ ∇h(x, t)

]2 + η(x, t) (1)

Here h(x, t) and Q(x, t) are the height of the landscape and the flow rate at position x and time t, η(x, t) is a noise term and α is a parameter which will be related to γ (see below). The constant term u on the right-hand-side of Eq. (1) mimics the so-called geological uplift which is known to originate by tectonic forces. The second term represents a local diffusion term of strength ν, condensing both smoothing (rounding of hilltops) and sedimentation (filling of valleys) processes. The third term is non-local and corresponds to an erosion driven by the waterflow on a surface with Q(x, t) representing the flow through site x. The exponent α is a parameter of the model which was assumed to be equal to 1 in Ref. [5] and Ref. [7]. The last noise term is added to account for small-scale stochastic processes such as rainfall fluctuations. In the attempt of extracting the basic features associated with the third term, one can ignore both the diffusion and noise terms (the accuracy of this approximation has been discussed in Ref. [5] ). In this case, a stationary state is obtained when the uplift term balances the flow dependent ero- sional term, that is when

|∇h(x, t)| ∼ Qγ−1(x, t) (2)

where, for later convenience, we define γ = 1−α/2. Most of previous work was performed with α = 1 (γ = 1/2) (there are few exceptions but with aims and methodolo- gies different from ours, see e.g. Ref. [1]). In this case Eq. (2) is known as slope-area law [1], which is a robust em- pirical law verified by real field observational data. We shall see later on how α is related to the dissipated en- ergy functional. Here we note that if we assume (as we shall do hereafter) uniform rainfall (and no ground wa- ter), then Q(x, t) ∼ a(x, t) where a(x, t) is the area of the basin draining into point x at time t. The basic result of this dynamical equation is that the non-linear term is able to account for the correct relationship between local slope and water flowrate.

III. CRITICAL EXPONENTS AND SCALING

LAWS

River networks are a remarkable example of systems governed by scaling laws. Although the concept of power and scaling laws has been known to the hydrologists for half a century, it was not until recently that this con- cept was put into a well defined framework [4], [18] in analogy with conventional non-equilibrium critical phe- nomena where power laws are associated with critical exponents. For the sake of completeness, we shall briefly review here a few of the central laws appearing in river networks which we regards as the most fundamental. Hack observed that the total area a draining into a

given point and the upstream length l going from that point to the source though the path of maximum water flow, were not independent but related by [21]

l ∼ ah (3)

where h is often referred as Hack’s exponent and it ranges in the interval [0.5 − 0.6] in real rivers. The distribu- tions of drainage areas and upstream lengths also follow a power law. Within the context of finite-size scaling, these may be written as [3,4]

p(a, L) = a−τf(a/Lφ) (4)

for areas and

pi(l, L) = l−ψg(l/Ldl) (5)

for lengths. In Eq. (4) p(a, L) is the distribution of drainage areas on a basin of size L and f(x) is a finite size function. It defines two critical exponents τ and φ which are not independent [4]. Similarly pi(l, L) is the distribu- tion of the upstream lengths on a basin of size L, g(x) is a finite size function and there exists a relation between ψ and dl. Note that φ and dl define the “typical” area (atyp ∼ L

φ) and length (ltyp ∼ L dl). Finally we remark

that one usually distinguishes between self-affine basins

2

(dl = 1, φ < 2) and self-similar basins (dl > 1, φ = 2) [22] and in both cases only one exponent out of four is independent. They are related via scaling relations as [4], [18], [22]

h = dl φ

τ = 2− h ψ = 1

h (6)

IV. ZERO-TEMPERATURE OCN: LOCAL

MINIMA

Within the context of the OCNs, the “dissipated en- ergy” of a river network at a given time, can be written, in continuum notations, as [14]:

E(t) =

∫ dx |∇h(x, t)|Q(x, t) (7)

The local gradient |∇h(x, t)| and the water flow rate Q(x, t) are expected to satisfy Eq. (2). Hence one gets:

E(t) =

∫ dx Qγ(x, t) (8)

We note that although the energy expression (8) has been derived in the context of OCN [2], in principle it could be defined in the evolution equation (1) as well, thus providing a way of monitoring the rate of approach to steady state. The basic assumption of the OCN is that there is a

tendency of any real river basin, to assume a configura- tion which minimizes Eq. (8). A natural question arising is the characterization of the local and absolute minima of E. It was shown [17,18], that the absolute minimum of E is insensitive to a variation of the value of γ in the interval 1/2 ≤ γ < 1 [20] and then h = 1/2, τ = 3/2 and dl = 1. The analogous issue in the case of local minima is much less clear. In fact, when γ = 0.5, numerical work [18] suggests that there exists a set of local minima which presumably corresponds to a new universality class that is the relevant one for real rivers. Indeed only exceptional events are able to radically modify river courses and so local minima of dissipated energy trap the system with high probability and therefore dominate the statistics. It is then a vital issue to discriminate whether or not

local minima configurations are indeed related to a well defined and robust universality class independent on γ.

V. RESULTS

In this section, we describe numerical procedures and results in detail for each case. All calculations are car- ried out on a L × L square lattice with periodic bound- ary conditions in one direction which we identify as the transverse direction. Multiple outlets are allowed in the

outflowing longitudinal direction (West side in the fig- ures) whereas an infinite wall is set up on the opposite side (East side). All averages considered in the statis- tics are carried out only over the river with the largest flow. It is worth stressing that the choice of considering only the maximum river in a multiple outlets environ- ment corresponds to consider the statistical behavior of rivers that are in competition to drain a given region. On the other hand, a single river within a given region is more appropriate if geological constraints are known to exist. Typically 8 nearest-neighbors (nn) - 4 associated with

the square lattice and 4 associated with the two diago- nal directions - are allowed. A somewhat more restricted choice considers only the 4 natural nn associated with the square lattice structure. In view of universality one would expect that the details of the lattice structure should not matter. This second choice has the considerable advan- tage of being less time consuming for numerical purposes. For this reason, although all following figures display pat- terns obtained with 8 nn, the results reported in this work are obtained with statistics based on 4 nn. In one example, we have explicitly checked that the outcomes using the two choices are consistent within the statistical errors. Finally, in all our networks, the drainage area a(x, t) is

computed, at each time step, in a standard way according to [1]:

a(x, t) = ∑ y(x)

a(y, t) + 1 (9)

where y(x) denotes all sites y which drain into x and the last term on the right hand side represents a unit rainfall input on each site at all times.

A. Initial condition

There are many possible initial conditions that can be used in numerical analysis of river networks. A popu- lar one among hydrologists is a deterministic comb-like structure [1]. However this initial condition suffers of various drawbacks [23] and its use in the multiple outlets case, such the one treated here, appears to be somewhat inconvenient. Hence we consider here two other phys- ically reasonable choices, namely a spanning tree (ST) and a Scheidegger network (SN). Although STs are well known in statistical physics mainly due to their relation with the q → 0 limit of the Potts model [24], only recently their topological properties have been studied in details [25]. A suitable variation of spanning trees has even been proposed as a topological model for river networks [22], [26]. We have generated STs with multiple outlets by us- ing an adapted Broder’s algorithm (see Ref. [25] and ref- erences therein). We have considered sizes ranging from L = 32 to L = 512 and statistics based on a number of configurations ranging from 500 to 100 respectively. We

3

have then performed a finite size analysis described in the next subsection to extract the most reliable values of the exponents. Our best result for the exponents are τ = 1.378± 0.002, ψ = 1.596± 0.003, h = 0.633± 0.003, dl = 1.25± 0.01, φ = 2.00 ± 0.01 in excellent agreement with the exact results τ = 1.375, ψ = 1.6, h = 0.625, dl = 1.25 [25]. This also provides a good test on the quality of our data analysis. Note that the exact value of φ, albeit not known, can be derived from the knowl- edge of h using Eq. (6). This predicts φ = 2.0 in perfect agreement with our computed value. We also note that the obtained results are nearly identical to the one pre- dicted and numerically generated for a single outlet [26]. A ST is a self-similar network in that it is undirected

and isotropic. A quite different choice (directed and anisotropic and hence self-affine) is a SN. This was again proposed as a topological model for river networks (see e.g. [1]) and the τ exponent has been exactly determined via a mapping to a one-dimensional model for mass ag- gregation [27]. A similar mapping to a diffusion-reaction model also provides the solution in general dimension- ality [28]. Our best results for the critical exponents with the same statistics as above are τ = 1.337± 0.003, ψ = 1.49 ± 0.02, h = 0.69 ± 0.02, dl = 0.96 ± 0.01, φ = 1.53 ± 0.02, which are in excellent agreement with the expected values (τ being exact, the others determined via Eq. (6)), that is τ = 4/3, ψ = 1.5, h = 0.75, dl = 1, φ = 1.5. In Figs.1 and 2, we depict typical patterns for a ST

and a SN respectively. The presence of multiple outlets magnifies the main difference between ST and SN. A ST is constructed in such a way that total freedom is given to the meandering of streams, the only constraint being that they all terminate on the same line (West side in the Figures). Typically this allows the formation of a main big river of size considerably larger with respect to the others. For SNs, the East-West preferred direction prevents the forming of such big river and many smaller rivers are usually present as shown in Fig.2.

B. The Dynamical Equation: A Self-Consistent

solution

As we discussed in Sec. II, we seek the stationary states of the following simplified equation:

∂h(x, t)

∂t = u−Qα(x, t)

[ ∇h(x, t)

]2 + η(x, t) (10)

where α = 2(1 − γ). The stationary averaged states of Eq. (10) are expected to conform to Eq. (2). Despite the apparent simplicity of the equation, an explicit numerical solution proves to be rather slow [5]. The reason for this can be traced back to the particular form of the erosion term. According to Eq. (10) only sites with a non-

negligible combination Qα(x, t) [ ∇h(x, t)

]2 will affect the

change of the pattern. As it turns out, this yields a long

time transient during which the elevation profile evolves very slowly. Since we are mainly interested in stationary states,

there is a way out from this situation [5]. The main idea is to start with an arbitrary network (e.g. a ST or a SN) and recursively construct the heights starting from the outlets with the aid of Eq. (2). From the derived landscape, a new network (in general different from the original one) can then be obtained by assuming that at each site the outward direction is along the steepest de- scent path and using Eq.(9). The noise term in Eq. (10) is mimicked by the unity term in Eq. (9). The proce- dure can then be iterated until self-consistency is finally achieved. The final configuration is, by definition, a sta- tionary state of Eq. (10). The convergence of the above procedure as a function

of the number of iterations is reported in Fig. 3 for various values of γ. Two remarks are in order. First we note that the dissipated energy as obtained from the above SC procedure does not have any physical meaning in the transient state. In other words the definition of “time” for the independent variable appearing in Fig. 3 should be considered only as a short-hand notation for “number of iterations”. Second it is apparent a the value γ = 0.5 seems to be the one with the slowest convergence ratio. This is probably due to the particular role played by the value γ = 0.5 as it will become clear shortly. Figs.4-8 depict typical patterns obtained on changing

the parameter γ in the interval [0, 1]. The effect of the value of γ on the aggregation pattern is evident. As γ increases from 0 to 1, single big rivers draining the entire basin in a snake-like form are less and less favored. The pattern for γ = 0 has a strong memory of the original initial tree. The overall effect of the SC procedure is to disfavor long meandering of the streams thus providing a self-affine character to the final tree. The main river be- comes rather straight for γ = 0.25 and the whole pattern appears to be more symmetric with respect to the case of γ = 0. A noteworthy feature is that this tendency is inverted as γ → 0.5 and returns back for γ = 0.75. As γ → 1 rivers become very directed as one would expect on the ground that this limit corresponds to the Schei- degger network [18]. The fact that γ = 0.5 most closely resembles real rivers is a reflection of the natural selec- tion by the erosional processes of this value of γ in terms of Eq. (2) as it has already been well documented in the literature (see e.g. [1]). In our numerical estimates of the exponents, for sim-

plicity, we shall restrict our attention to the half region [0, 0.5] . We also note that this is the region inaccessible to the analytical scheme of Ref. [17]. For a more accurate evaluation of the critical expo-

nents τ and ψ it proves convenient to introduce the inte- grated probabilities:

P (a, L) =

∫ a 0

da′p(a′, L) = a1−τF ( a

Lφ ) (11)

and

4

Π(l, L) =

∫ l 0

dl′pi(l′, L) = l1−ψG( l

Ldl ) (12)

where F (x) and G(x) are related to f(x) and g(x) de- fined in Eqs. (4) and (5), in an obvious way. An efficient way of computing the exponents is through the so-called “effective” (sometime also referred to as “running”) ex- ponents. In the present case they are defined as:

τ(a) = 1− ∂ logP (a, L)

∂ a (13)

and

ψ(l) = 1− ∂ logΠ(l, L)

∂ l . (14)

One then obtains an effective exponent for each value of the independent variable (a or l). Fig. 9 shows one typical result on a 256× 256 lattice.

We can divide the obtained values roughly in four regions. The first region (1 ≤ a < 10) corresponds to the region of no scaling. Small rivers belonging to the second region (10 ≤ a < 100) have an exponent close to the absolute minimum value τ = 1.5. This is consistent with the pic- ture of typical rivers (see Fig. 6) where small rivers dis- play a marked straightness similar to the one of the abso- lute minimum [18] and it means that they quickly assume configurations consistent with the absolute minimum (in- dependent of the initial conditions). Larger areas are as- sociated with larger rivers which have longer memory of the initial condition (ST in the present case). Hence the corresponding exponent is sensibly smaller (closer to the ST value 1.38). The last region corresponds to the finite size cut-off and must be discarded. After discarding the first and forth regions the ob-

tained values can then be grouped into local bins and a local average exponent can be associated with each of them. Statistical fluctuations within each box then yield an estimate of error bars. This provides our best esti- mate of the exponent for each value of L and a simple 1/L extrapolation is then carried out to extract the final values. This is depicted in Fig. 10 and the corresponding best estimates of this method are reported in Table I. An analogous procedure leads to the best estimates for ψ as reported again in Table I. Sizes and statistics are identical to those considered for the initial conditions and hence these simulations are rather time consuming. One can notice a weak dependence on γ for both τ and

ψ. One way of computing φ and dl is through the collapse plots of the probabilities P (a, L) and Π(l, L). However we have found that a satisfactory collapse can achieve only within a limited range of the appropriate variable (a/Lφ and l/Ldl). Hence we have opted for an alternative scheme hinging on the calculation of the following ratios:

M qa (L) ≡

⟨ aq+1

⟩ a

〈aq〉a ∼ Lφ q = 1, 2, ... (15)

where averages are over the probability densities p(a, L) (of the maximum river) and the L dependence is straight- forward [4]. A similar relationship holds for the lengths:

M ql (L) ≡

⟨ lq+1

⟩ l

〈lq〉l ∼ Ldl q = 1, 2, ... (16)

The results for q = 1, 2, 3 are reported in Table III for the exponent φ and in Table IV for the exponent dl. Sur- prisingly the values for φ are consistently larger that the space-filling value φ = 2 which is expected to be the up- per bound for this exponent. This is probably due to a deficiency of this procedure during a cross-over from a self-similar regime (the initial ST) to a self-affine pattern (the final configuration). Indeed we shall see later on that this feature is not present when one starts with a self-affine network (e.g. a SN) from the outset. A seem- ingly large value of dl is appearing probably due to the same reason. Overall, one can notice a rather weak (if any) dependence on γ and almost no dependence on q. Finally Hack’s exponent h was computed from its def-

inition Eq. (3) using an effective exponent method to be described below. Let us assume a power-law dependence for a generic

function as given by

f(x) ∼ xθ θ > 0 (17)

On integrating f(x) (between lower and upper limits, say x0 and x), we find

θ =

[ xf(x) − x0f(x0)

] ∫ x x0 dzf(z)

− 1. (18)

The effective exponent obtained using Eq.(18) in Eq.(3) with f(x) ≡ a, x ≡ l, and θ ≡ 1/h, can then be ana- lyzed with the same procedure (local average plus 1/L extrapolation) outlined before. We note that for a we have used an averaged value over all areas correspond- ing to the same length. The final results are reported in Table I and one can see that there is an overall good agreement with scaling laws.

C. OCN

The minimization of the energy functional Eq. (8) goes through an algorithm akin to the one exploited by Sun et. al [16]. It is based on the following steps:

1) An initial configuration (ST or SN) is generated and its dissipated energy is computed according to Eq. (8). By definition, this initial network has no loops.

2) A link is randomly selected and its local outflow is also randomly chosen.

5

3) This new configuration is tested for loop creation. If a loop has been created, the configuration is rejected and step 2) is repeated.

4) The energy of the new candidate configuration is com- puted. If it is smaller than the previous one, the new configuration is accepted, otherwise it is re- jected.

5) Steps 2)-4) are repeated until the energy does not change within a given tolerance.

The final configuration is regarded as a local minimum. This scheme is patterned after a standard Metropolis al- gorithm at zero temperature, and averages are over many different configurations (ranging from 500 at L = 32 to 100 at L = 256). Fig. 11 depicts the dissipated energy per unit of length as a function of the convergence “time” (i.e. the number of total iteration of the algorithm). The similarity with patterns obtained from the SC procedure is evident. Here, however, the typical convergence time is much longer than before. Critical exponents are computed with the same pre-

scription given in the previous subsection. In Fig. 12 we show the resulting 1/L extrapolation. The final best estimates are reported in Table I. Once more, a weak dependence on γ (τ increases as γ

increases) can be noticed. All other exponents are consis- tent with scaling relations Eq. (6). It is worth stressing that exponents are nearly consistent (at the edge of sta- tistical errors) with those previously obtained from the dynamical equation. Regarding the exponents φ and dl, they can be found in

Tables III and IV respectively. Here too the same feature, discussed in VB in connection with the exponent values of φ and dl, applies.

D. Independence of the initial condition

Our final task is a test of the sensitivity of critical ex- ponents to the initial conditions. To this aim we have changed the initial condition from a ST to a SN for both the dynamical equation and the OCN. The main differ- ence between the two initial conditions is that while the latter is a directed network, the former is not. Table II reports the comparison for the case γ = 0.5, and Fig. 13 depicts the network resulting from the SC scheme, for a typical final state. Despite the obvious memory of the initial network (see Fig. 2), the typical distance among big rivers is larger than the initial one. This in turn yields a higher value for H and hence non-trivial exponents. It is remarkable that although these patterns appear

to be considerably different from those obtained when starting with a ST, the two sets of critical exponents are nevertheless consistent within the statistical errors. As a final comment, we find in this case values of φ

and dl in very good agreement with scaling predictions. The results for different ratios q in the case γ = 0.5 for

both the SC and the OCN, along with the comparison with the corresponding values stemming from STs are reported in Tables V and VI respectively.

VI. CONCLUSIONS

In this paper we have addressed the issue of the exis- tence and robustness of the universality class associated with landscapes corresponding to local minimal energies. To this aim we have first extended previous studies

for both the SC solution of a Langevin equation and the OCNs variational methods. Higher sizes and statistics have been exploited in both cases and, (to our knowledge) for the first time, the most physical procedure of basing the statistics only on the river with largest flow, has been implemented. Second, we have monitored the dependence of critical

exponents on a parameter γ associated with the slope- area law in one case (SC) and with the dissipated energy in the other (OCN). Our results give compatible critical exponents between SC and OCN within the error bars, but a weak and similar dependence on γ appears in both models. Finally we have tested the stability of the obtained

results for both the SC and the OCN with respect to changes in the initial conditions. Although the obtained final patterns display a dependence on initial conditions, critical exponents appear to be insensitive to this depen- dence. As a by-product of our investigation, we have found

that the SC solution of the dynamical equation is a very powerful method to investigate river networks as it is ca- pable of providing useful informations on the stationary state in a simple and physical way. Another interesting point, from a numerical point of view, is that this pro- cedure typically achieves convergence much faster with respect to OCN scheme. In view of our results we can now summarize the argu-

ments favoring and disfavoring the appearance of a new universality class associated with local minima. As we mentioned in our discussion of Fig. 10 the typical evolu- tion of a network appears to depend on the considered length scale. Small rivers very quickly settle to their final state whereas much longer time is required to large rivers to forget their initial conditions. This is also reflected by the difficulty in collapsing the distribution probabilities of both a and l into a single plot for a reasonably ex- tended range of the corresponding variable. It is then possible that although only two universality classes (one associated with the initial condition and the other to the absolute minimum) are present, an intermediate univer- sality class sets on due to both the averaging over dif- ferent regimes and the difficulty of reaching the absolute minimum. On the other hand, in support of the existence of a new

universality class, we cite the fact that critical exponents

6

are different from those of both ST and SN – on start- ing with these initial conditions, critical exponents in the final configurations clearly deviate from their initial val- ues. Furthermore all exponents are found to be robust and obey to scaling relations summarized in Sec.III. Overall we believe that the evidence suggested by our

results favors this second possibility which has also been hinted into a different context [29]. In this respect the weak γ dependence of critical exponents remains unex- plained. It would be interesting to generalize the present cal-

culation in two aspects. For the dynamical equation, it would be instructive to tackle the problem of the ex- plicit solution of Eq. (1). This route has the advantage that Eq. (2) (the key relation between erosional process and network topology) can be then derived rather than assumed as in the self-consistent procedure. Similarly, a parallel calculation could also be implemented in the OCN framework, upon starting with the more general expression given in Eq. (7).

ACKNOWLEDGMENTS

It is a pleasure to thank Jayanth Banavar, Amos Maritan and Andrea Rinaldo for many enlighting dis- cussions, innumerable suggestions and a critical reading of the manuscript. Financial support from the Italian MURST (Ministero dell’Universita´ e della Ricerca Sci- entifica e Tecnologica) and INFM (Istituto Nazionale di Fisica della Materia) is gratefully acknowledged.

FIG. 1. A spanning tree on a 64× 64 square lattice with 8 nearest-neighbors, 4 in the Nord-South and East-West direc- tion and 4 along the diagonal directions. The outlet is on the West side (outflowing), on the East side there is an infinite wall, whereas there are periodic boundary conditions in the North-South direction. The thickness of the line at each point is proportional to the flow through that point. The seeming loops are just an artifact of the drawing.

FIG. 2. A Scheidegger network with L = 64. Boundary conditions are the same as in Fig.1. The directness of the network provides a privileged East-West direction.

7

0 500 1000 t

2500

3500

4500 E

/L

γ=0.00 γ=0.25 γ=0.50

FIG. 3. Dissipated energy per unit of length, as a function of the number of iterations (time t) of the self-consistent solu- tion of the dynamical equation, on starting with the spanning tree of Fig 1. The parameter γ is the value appearing in Eq. (2). The system size is L = 512 and each point of the curve is an average over all configurations which have gone at least that far in the number of iterations. All quantities reported in this and following figures are dimensionless.

FIG. 4. A typical network obtained as a final output of the self-consistent procedure described in Sec. VB. Here L = 64, γ = 0, and the initial condition is the spanning tree of Fig. 1.

FIG. 5. Same as in Fig. 4, with γ = 0.25.

FIG. 6. Same as in Fig. 4, with γ = 0.5.

8

FIG. 7. Same as in Fig. 4, with γ = 0.75.

FIG. 8. Same as in Fig. 4, with γ = 1.00.

0 200 400 600 800 1000 a

1.0

1.1

1.2

1.3

1.4

1.5

1.6

1.7

1.8

1.9

2.0

τ

FIG. 9. Effective exponent τ as a function of the area a in the stationary state of the dynamical equation with γ = 0.5. Here L = 256 and the initial condition is the spanning tree of Fig. 1.

0.00 0.01 0.02 0.03 0.04 1/L

1.0

1.1

1.2

1.3

1.4

1.5

1.6

1.7

1.8

1.9

2.0

τ

γ=0.00 γ=0.25 γ=0.50

FIG. 10. Finite size 1/L extrapolation for the value of τ obtained with the dynamical equation for the cases γ = 0, 0.25, 0.5. In all cases the initial condition is the span- ning tree of Fig. 1.

9

0 50000 100000 t

500

1000

1500

2000

2500 E

/L

γ=0.25 γ=0.50

FIG. 11. Dissipated energy per unit of length, as a function of the number of iterations in the OCN procedure. Here the size is L = 256 and again each point of the curves are averaged over all configurations which have reached that time t. The initial condition is the spanning tree of Fig. 1.

0.00 0.01 0.02 0.03 0.04 1/L

1.0

1.1

1.2

1.3

1.4

1.5

1.6

1.7

1.8

1.9

2.0

τ

γ=0.25 γ=0.50

FIG. 12. Finite size 1/L extrapolation for the value of τ obtained with the OCN (T = 0) for the cases γ = 0.25, 0.5. In both cases the initial condition is the spanning tree of Fig. 1.

FIG. 13. A typical 64 × 64 network obtained as a final output of the self-consistent procedure described in Sec. VB upon starting with a the Scheidegger network of Fig. 2.

TABLE I. Critical exponents τ , ψ and h as a function of γ for both the Self-Consistent (SC) procedure and the optimal channel networks (OCN) at zero temperature. The values in the parenthesis are those obtained from the computed value of τ and using the scaling relations. For γ = 0, the OCN exponents reported are the exact ones corresponding to the initial conditions (a spanning tree in the present case) since the energy is a constant.

γ τOCN τSC ψOCN ψSC hO 0.00 1.38 1.40± 0.01 1.6 1.69± 0.03 (1.67) 0.625 0.25 1.42± 0.01 1.42± 0.01 1.71± 0.07 (1.72) 1.68± 0.05 (1.72) 0.64± 0 0.50 1.44± 0.01 1.46± 0.01 1.8± 0.1 (1.79) 1.82± 0.05 (1.85) 0.61± 0

TABLE II. Critical exponents τ , ψ and h for γ = 0.5 as a function of the initial conditions. As in the text, ST and SN stand for Spanning Trees and Scheidegger Network respec- tively. SC and OCN have the same meaning as above. Values in parenthesis are scaling predictions.

τOCN τSC ψOCN ψSC hOCN ST 1.44 ± 0.01 1.46 ± 0.01 1.8 ± 0.1 (1.79) 1.82± 0.05 (1.85) 0.61± 0.03 SN 1.44 ± 0.03 1.43 ± 0.02 1.8 ± 0.2 (1.79) 1.7± 0.1 (1.75) 0.61± 0.03

10

TABLE III. Critical exponent φ stemming from both the SC procedure and the OCN scheme, as a function of the value of γ and of the value of the order q of ratios defined in Eqs. (13) and (14).

q γ = 0.00 (SC ) γ = 0.25 (SC ) γ = 0.50 (SC ) γ = 0.00 (OCN) γ = 0.25 (OCN) γ = 0.50 (OCN)

1 2.1 ± 0.1 2.2± 0.1 2.2 ± 0.1 - 2.0± 0.1 2.1± 0.1 2 2.1 ± 0.1 2.2± 0.1 2.2 ± 0.1 - 2.0± 0.1 2.1± 0.1 3 2.1 ± 0.1 2.2± 0.1 2.3 ± 0.1 - 2.0± 0.1 2.0± 0.1

TABLE IV. Critical exponent dl obtained from both the SC procedure and the OCN scheme, as a function of the value of γ and of the value of the order q of ratios (13) and (14) .

q γ = 0.00 (SC ) γ = 0.25 (SC ) γ = 0.50 (SC ) γ = 0.00 (OCN) γ = 0.25 (OCN) γ = 0.50 (OCN)

1 1.3 ± 0.1 1.3± 0.1 1.1 ± 0.1 - 1.2± 0.1 1.3± 0.1 2 1.3 ± 0.1 1.4± 0.1 1.2 ± 0.1 - 1.3± 0.1 1.3± 0.1 3 1.3 ± 0.1 1.4± 0.1 1.2 ± 0.1 - 1.3± 0.1 1.3± 0.1

TABLE V. Comparison between exponents φ as obtained from both the SC procedure and the OCN scheme, as a func- tion of the value of the initial conditions and of the value of the order q of ratios (13) and (14). As in the text, ST stands for Spanning Trees and SN for Scheidegger Network. Here γ = 0.5.

q φST (SC ) φSN (SC ) φST (OCN) φSN (OCN)

1 2.2 ± 0.1 1.8± 0.1 2.1± 0.1 1.7± 0.1 2 2.2 ± 0.1 1.8± 0.1 2.1± 0.1 1.7± 0.1 3 2.3 ± 0.1 1.8± 0.1 2.0± 0.1 1.7± 0.1

TABLE VI. Comparison between exponents dl resulting from both the SC procedure and the OCN scheme, as a func- tion of the value of the initial conditions and of the value of the order q of ratios (13) and (14). As in the text, ST stands for Spanning Trees and SN for Scheidegger Network. Here γ = 0.5.

q dlST (SC ) dlSN (SC ) dlST (OCN) dlSN (OCN)

1 1.1 ± 0.1 0.9± 0.1 1.3 ± 0.1 0.9± 0.1 2 1.2 ± 0.1 1.0± 0.1 1.3 ± 0.1 1.0± 0.1 3 1.2 ± 0.1 1.0± 0.1 1.3 ± 0.1 1.0± 0.1

[1] I. Rodrigues-Iturbe and A. Rinaldo, Fractal River Basins: Chance and Self-Organization (Cambridge Uni- versity Press, Great Britain, 1997).

[2] A. Rinaldo, I. Rodrigues-Iturbe, R. Rigon and R.L. Bras, Phys. Rev. Lett. 70, 822 (1993).

[3] P. Meakin, J. Feder, and T. Jossang, Physica A 176, 409 (1991).

[4] A. Maritan, A. Rinaldo, R. Rigon, A. Giacometti and I. Rodrigues-Iturbe, Phys. Rev. E 53, 1510 (1996).

[5] J.R. Banavar, F. Colaiori, A. Flammini, A. Giacometti, A. Maritan and A. Rinaldo, Phys. Rev. Lett. 78, 4522 (1997).

[6] K. Sinclair and R. C. Ball, Phys. Rev. Lett. 76, 3360 (1996).

[7] E. Somafai and L. M. Sander, Phys. Rev. E 56, R5 (1997).

[8] P. S. Dodds and D. H. Rothman, Phys. Rev. E 59, 4865 (1999).

[9] L.B.Leopold and W.B.Langbein, U.S. Geol. Surv. Prof. bf 500-A,1 (1962).

[10] M. Marsili, A. Maritan, F. Toigo and J. R. Banavar, Rev. Mod. Phys. 68, 963 (1996).

[11] R. Pastor-Satorras and D. H. Rothman, Phys. Rev. Lett. 19, 4349 (1998).

[12] F. Me´tiver, Phys. Rev. E 60, 5827 (1999) [13] A. Giacometti, A. Maritan, and J. R. Banavar, Phys.

Rev. Lett. 75, 577 (1995). [14] I. Rodriguez-Iturbe, A. Rinaldo, R. Rigon, R. L. Bras,

and E. Ijjasz-Vasquez, Water Resour. Res. 28, 1095 (1992); I. Rodriguez-Iturbe, A. Rinaldo, R. Rigon, R. L. Bras, and E. Ijjasz-Vasquez, Geophys. Res. Lett. 19, 889 (1992); A. Rinaldo et. al, Water Resour. Res. 28, 2183 (1992).

[15] S. Kirkpatrick, C.D. Gelatt, and M. P. Vecchi, Science 220, 671 (1983).

[16] T. Sun, P. Meakin and T. Jossang, Phys. Rev. E 49, 4865 (1994).

[17] A. Maritan, F. Colaiori, A. Flammini, M. Cieplak and J. R. Banavar, Science 272, 984 (1996).

[18] F. Colaiori, A. Flammini, A. Maritan and J. R. Banavar, Phys. Rev. E 55, 1298 (1997).

[19] The Peano basin is discussed e.g. in B.B. Mandelbrot The Fractal Geometry of Nature, (Freeman, New York 1983). Some topological properties in the river network framework have been addressed in A. Marani, R. Rigon and A. Rinaldo, Water Resour. Res. 27, 3041 (1991).

[20] We note that γ = 0 and γ = 1 correspond to known and exactly solvable models (spanning trees and Scheidegger network respectively) and that the only physically ac- ceptable values of γ lie between these two limiting cases.

[21] J. T. Hack, U.S. Geological Survey Professional Papers 294-B, 45 (1957).

[22] M. Cieplak, A. Giacometti, A. Maritan, A. Rinaldo, I. Rodrigues-Iturbe and J. R. Banavar, J. Stat. Phys. 91, 1 (1998).

11

[23] G. Caldarelli, A. Giacometti, A. Maritan, I.Rodrigues- Iturbe and A. Rinaldo, Phys. Rev. E 55, R4865 (1997).

[24] F. Y. Wu, Rev. Mod. Phys. 54, 235 (1982). [25] S. S. Manna, D. Dhar and S. N. Majumadar, Phys. Rev.

A 46, R4471 (1992). [26] S.S. Manna and B. Subramanian, Phys. Rev. Lett. 76,

3460 (1996). [27] M. Takayasu, H. Takayasu and G. Huber, J. Stat. Phys.

65, 725 (1991). [28] M. R. Swift, F. Colaiori, A. Flammini, A. Maritan, A.

Giacometti and J. R. Banavar, Phys. Rev. Lett. 79, 3278 (1997).

[29] B. Tadic, Phys. Rev. E 58, 168 (1998).

12

iv :c

on d-

m at

/0 00

84 24

v1 [

co nd

-m at.

sta t-m

ec h]

2 9 A

ug 20

00 Local minimal energy landscapes in river networks

Achille Giacometti INFM Unita´ di Venezia, Dipartimento di Scienze Ambientali, Universita´ di Venezia

Calle Larga Santa Marta DD2137, I-30123 Venezia-Italy

(February 1, 2008)

The existence and stability of the universality class as- sociated to local minimal energy landscapes is investigated. Using extensive numerical simulations, we first study the de- pendence on a parameter γ of a partial differential equation which was proposed to describe the evolution of a rugged land- scape toward a local minimum of the dissipated energy. We then compare the results with those obtained by an evolution scheme based on a variational principle (the optimal channel networks). It is found that both models yield qualitatively similar river patterns and similar dependence on γ. The ag- gregation mechanism is however strongly dependent on the value of γ. A careful analysis suggests that scaling behaviors may weakly depend both on γ and on initial condition, but in all cases it is within observational data predictions. Con- sequences of our results are finally discussed and the most plausible scenario is presented.

92.40.Fb,64.60.Ht,05.60.+w

I. INTRODUCTION

Understanding the development of a landscape in terms of fundamental mechanical principles is a formidable task [1]. In spite of this high complexity, re- cent theoretical studies have resulted in considerable pro- gresses by considering the issue from a viewpoint anal- ogous to the one taken in conventional critical phenom- ena [2–8] where simple models are exploited to identify universality classes. The main idea of this approach is indeed to focus on few fundamental ingredients which, in the spirit of critical phenomena, are expected to provide a reasonable description of the large-scales properties of the network. The remarkable properties of river networks have been

known for some time and they can be condensed in few phenomenological scaling laws which have been con- firmed to hold in the observational data [9]. While these laws do not explain the underlying physical mech- anisms, they nevertheless provide guidelines for their search. Hence any newly proposed model for river net- works ought to be tested against these laws. In the lan- guage of critical phenomena, those scaling laws can be used to derive critical exponents and thus discriminate among different universality classes. Among such laws a fundamental role in the physics of the erosion of a land- scape due to the flow of water over it, is played by the slope-area law (see Ref. [1] for a review).

A few Langevin-like equations have been recently proposed for describing the evolution of the landscape under the effect of erosional processes. In Ref. [5], reparametrization invariance arguments [10] were used to derive a dynamical equation which yields the slope- area law as a stationary state. The same equation was obtained by Somafai and Sander [7] using Landau ar- guments. Other proposals have also been advanced [6], [11–13]. The equation proposed in Ref. [5], was studied both

analytically in 1+1 dimensions (one spatial and one tem- poral) and numerically using a self-consistent (SC) solu- tion in 2+1 dimensions, and it was shown to predict a fairly reasonable stationary state quite different from the starting network acting as initial condition. This type of analysis hinges on the observation that the network rep- resenting the river has an evolution mirroring that of the evolution profile. Other methods have also been devised which directly address the topological properties of the network itself, the best known being the optimal channel networks (OCN) [14]. This is a lattice model where a functional describing the dissipated energy is minimized in order to find the optimal configuration and it is based on the idea that, presumably, the erosional process tak- ing place on a landscape is driven by a strive for opti- mality. A simulating annealing procedure [15] has been implemented to this aim both at zero [2] and finite tem- perature [16]. While the latter is aimed to find the ab- solute minimum, the former is expected to display local minima only. Using exact bounds and finite-size scaling, Maritan et al ( [17,18]) showed that the absolute mini- mum belongs to the mean-field universality class, that, in turn, means that the corresponding network has a highly symmetric pattern with small rivers draining into bigger rivers in a predictable way (this network is akin to the Peano basin [19]). However this minimum is not eas- ily reachable in the space of all configurations, and one is then led to suspect that real rivers are better described by configurations related to local minima. We further note that stationary states of the aforementioned dynamical equations, are also expected to be associated with local minima for the same reason. In view of the above discussion, it is apparent that a

more complete description of the stability and the scaling behavior associated with these local minima would be de- sirable. The present work is an attempt in this direction. Our aims are threefold. Our first goal is the full char- acterization of the possible universality class associated with local minima. Using the SC numerical procedure

1

envisaged in Ref. [5], we extend previous numerical re- sults both in size and statistics and we use a finite size procedure for a more accurate estimate of the exponents. We then compare these exponents with those calculated using a zero-temperature simulation in the OCN frame- work with similar sizes and statistics. Our second goal is the generalization of the dynamical equation to include a tunable parameter which turns out to be related to the 0 ≤ γ ≤ 1 parameter considered in OCNs and which is responsible of the aggregation mechanism. In Ref. [18] it was shown that the absolute minimum is insensitive to a variation of γ in the interval 1/2 ≤ γ < 1 [20]. We find that although the final patterns display marked differences as a function of γ, critical exponents show a much smaller discrepancy, which our results indicate to be marginal (at the edge of error bars), both in the SC and OCN case. Finally we test the stability of our results against a variation in the initial conditions. We find that even in this case a marginal difference appears in the effective critical exponents. The paper is organized as follows. In Sec.II, we re-

view the dynamical equation whereas in Sec.III scaling laws and critical exponents are briefly recalled. Sec.IV contains the definition of OCN and our results are pre- sented in Sec.V. Finally Sec.VI contains some concluding remarks and future perspectives.

II. THE DYNAMICAL EQUATION

A rather general dynamical equation consistent with general principles, and capturing the physics of erosional processes occurring in river basin, is [5], [7]:

∂h(x, t)

∂t = u+ ν∇2h(x, t)−Qα(x, t)

[ ∇h(x, t)

]2 + η(x, t) (1)

Here h(x, t) and Q(x, t) are the height of the landscape and the flow rate at position x and time t, η(x, t) is a noise term and α is a parameter which will be related to γ (see below). The constant term u on the right-hand-side of Eq. (1) mimics the so-called geological uplift which is known to originate by tectonic forces. The second term represents a local diffusion term of strength ν, condensing both smoothing (rounding of hilltops) and sedimentation (filling of valleys) processes. The third term is non-local and corresponds to an erosion driven by the waterflow on a surface with Q(x, t) representing the flow through site x. The exponent α is a parameter of the model which was assumed to be equal to 1 in Ref. [5] and Ref. [7]. The last noise term is added to account for small-scale stochastic processes such as rainfall fluctuations. In the attempt of extracting the basic features associated with the third term, one can ignore both the diffusion and noise terms (the accuracy of this approximation has been discussed in Ref. [5] ). In this case, a stationary state is obtained when the uplift term balances the flow dependent ero- sional term, that is when

|∇h(x, t)| ∼ Qγ−1(x, t) (2)

where, for later convenience, we define γ = 1−α/2. Most of previous work was performed with α = 1 (γ = 1/2) (there are few exceptions but with aims and methodolo- gies different from ours, see e.g. Ref. [1]). In this case Eq. (2) is known as slope-area law [1], which is a robust em- pirical law verified by real field observational data. We shall see later on how α is related to the dissipated en- ergy functional. Here we note that if we assume (as we shall do hereafter) uniform rainfall (and no ground wa- ter), then Q(x, t) ∼ a(x, t) where a(x, t) is the area of the basin draining into point x at time t. The basic result of this dynamical equation is that the non-linear term is able to account for the correct relationship between local slope and water flowrate.

III. CRITICAL EXPONENTS AND SCALING

LAWS

River networks are a remarkable example of systems governed by scaling laws. Although the concept of power and scaling laws has been known to the hydrologists for half a century, it was not until recently that this con- cept was put into a well defined framework [4], [18] in analogy with conventional non-equilibrium critical phe- nomena where power laws are associated with critical exponents. For the sake of completeness, we shall briefly review here a few of the central laws appearing in river networks which we regards as the most fundamental. Hack observed that the total area a draining into a

given point and the upstream length l going from that point to the source though the path of maximum water flow, were not independent but related by [21]

l ∼ ah (3)

where h is often referred as Hack’s exponent and it ranges in the interval [0.5 − 0.6] in real rivers. The distribu- tions of drainage areas and upstream lengths also follow a power law. Within the context of finite-size scaling, these may be written as [3,4]

p(a, L) = a−τf(a/Lφ) (4)

for areas and

pi(l, L) = l−ψg(l/Ldl) (5)

for lengths. In Eq. (4) p(a, L) is the distribution of drainage areas on a basin of size L and f(x) is a finite size function. It defines two critical exponents τ and φ which are not independent [4]. Similarly pi(l, L) is the distribu- tion of the upstream lengths on a basin of size L, g(x) is a finite size function and there exists a relation between ψ and dl. Note that φ and dl define the “typical” area (atyp ∼ L

φ) and length (ltyp ∼ L dl). Finally we remark

that one usually distinguishes between self-affine basins

2

(dl = 1, φ < 2) and self-similar basins (dl > 1, φ = 2) [22] and in both cases only one exponent out of four is independent. They are related via scaling relations as [4], [18], [22]

h = dl φ

τ = 2− h ψ = 1

h (6)

IV. ZERO-TEMPERATURE OCN: LOCAL

MINIMA

Within the context of the OCNs, the “dissipated en- ergy” of a river network at a given time, can be written, in continuum notations, as [14]:

E(t) =

∫ dx |∇h(x, t)|Q(x, t) (7)

The local gradient |∇h(x, t)| and the water flow rate Q(x, t) are expected to satisfy Eq. (2). Hence one gets:

E(t) =

∫ dx Qγ(x, t) (8)

We note that although the energy expression (8) has been derived in the context of OCN [2], in principle it could be defined in the evolution equation (1) as well, thus providing a way of monitoring the rate of approach to steady state. The basic assumption of the OCN is that there is a

tendency of any real river basin, to assume a configura- tion which minimizes Eq. (8). A natural question arising is the characterization of the local and absolute minima of E. It was shown [17,18], that the absolute minimum of E is insensitive to a variation of the value of γ in the interval 1/2 ≤ γ < 1 [20] and then h = 1/2, τ = 3/2 and dl = 1. The analogous issue in the case of local minima is much less clear. In fact, when γ = 0.5, numerical work [18] suggests that there exists a set of local minima which presumably corresponds to a new universality class that is the relevant one for real rivers. Indeed only exceptional events are able to radically modify river courses and so local minima of dissipated energy trap the system with high probability and therefore dominate the statistics. It is then a vital issue to discriminate whether or not

local minima configurations are indeed related to a well defined and robust universality class independent on γ.

V. RESULTS

In this section, we describe numerical procedures and results in detail for each case. All calculations are car- ried out on a L × L square lattice with periodic bound- ary conditions in one direction which we identify as the transverse direction. Multiple outlets are allowed in the

outflowing longitudinal direction (West side in the fig- ures) whereas an infinite wall is set up on the opposite side (East side). All averages considered in the statis- tics are carried out only over the river with the largest flow. It is worth stressing that the choice of considering only the maximum river in a multiple outlets environ- ment corresponds to consider the statistical behavior of rivers that are in competition to drain a given region. On the other hand, a single river within a given region is more appropriate if geological constraints are known to exist. Typically 8 nearest-neighbors (nn) - 4 associated with

the square lattice and 4 associated with the two diago- nal directions - are allowed. A somewhat more restricted choice considers only the 4 natural nn associated with the square lattice structure. In view of universality one would expect that the details of the lattice structure should not matter. This second choice has the considerable advan- tage of being less time consuming for numerical purposes. For this reason, although all following figures display pat- terns obtained with 8 nn, the results reported in this work are obtained with statistics based on 4 nn. In one example, we have explicitly checked that the outcomes using the two choices are consistent within the statistical errors. Finally, in all our networks, the drainage area a(x, t) is

computed, at each time step, in a standard way according to [1]:

a(x, t) = ∑ y(x)

a(y, t) + 1 (9)

where y(x) denotes all sites y which drain into x and the last term on the right hand side represents a unit rainfall input on each site at all times.

A. Initial condition

There are many possible initial conditions that can be used in numerical analysis of river networks. A popu- lar one among hydrologists is a deterministic comb-like structure [1]. However this initial condition suffers of various drawbacks [23] and its use in the multiple outlets case, such the one treated here, appears to be somewhat inconvenient. Hence we consider here two other phys- ically reasonable choices, namely a spanning tree (ST) and a Scheidegger network (SN). Although STs are well known in statistical physics mainly due to their relation with the q → 0 limit of the Potts model [24], only recently their topological properties have been studied in details [25]. A suitable variation of spanning trees has even been proposed as a topological model for river networks [22], [26]. We have generated STs with multiple outlets by us- ing an adapted Broder’s algorithm (see Ref. [25] and ref- erences therein). We have considered sizes ranging from L = 32 to L = 512 and statistics based on a number of configurations ranging from 500 to 100 respectively. We

3

have then performed a finite size analysis described in the next subsection to extract the most reliable values of the exponents. Our best result for the exponents are τ = 1.378± 0.002, ψ = 1.596± 0.003, h = 0.633± 0.003, dl = 1.25± 0.01, φ = 2.00 ± 0.01 in excellent agreement with the exact results τ = 1.375, ψ = 1.6, h = 0.625, dl = 1.25 [25]. This also provides a good test on the quality of our data analysis. Note that the exact value of φ, albeit not known, can be derived from the knowl- edge of h using Eq. (6). This predicts φ = 2.0 in perfect agreement with our computed value. We also note that the obtained results are nearly identical to the one pre- dicted and numerically generated for a single outlet [26]. A ST is a self-similar network in that it is undirected

and isotropic. A quite different choice (directed and anisotropic and hence self-affine) is a SN. This was again proposed as a topological model for river networks (see e.g. [1]) and the τ exponent has been exactly determined via a mapping to a one-dimensional model for mass ag- gregation [27]. A similar mapping to a diffusion-reaction model also provides the solution in general dimension- ality [28]. Our best results for the critical exponents with the same statistics as above are τ = 1.337± 0.003, ψ = 1.49 ± 0.02, h = 0.69 ± 0.02, dl = 0.96 ± 0.01, φ = 1.53 ± 0.02, which are in excellent agreement with the expected values (τ being exact, the others determined via Eq. (6)), that is τ = 4/3, ψ = 1.5, h = 0.75, dl = 1, φ = 1.5. In Figs.1 and 2, we depict typical patterns for a ST

and a SN respectively. The presence of multiple outlets magnifies the main difference between ST and SN. A ST is constructed in such a way that total freedom is given to the meandering of streams, the only constraint being that they all terminate on the same line (West side in the Figures). Typically this allows the formation of a main big river of size considerably larger with respect to the others. For SNs, the East-West preferred direction prevents the forming of such big river and many smaller rivers are usually present as shown in Fig.2.

B. The Dynamical Equation: A Self-Consistent

solution

As we discussed in Sec. II, we seek the stationary states of the following simplified equation:

∂h(x, t)

∂t = u−Qα(x, t)

[ ∇h(x, t)

]2 + η(x, t) (10)

where α = 2(1 − γ). The stationary averaged states of Eq. (10) are expected to conform to Eq. (2). Despite the apparent simplicity of the equation, an explicit numerical solution proves to be rather slow [5]. The reason for this can be traced back to the particular form of the erosion term. According to Eq. (10) only sites with a non-

negligible combination Qα(x, t) [ ∇h(x, t)

]2 will affect the

change of the pattern. As it turns out, this yields a long

time transient during which the elevation profile evolves very slowly. Since we are mainly interested in stationary states,

there is a way out from this situation [5]. The main idea is to start with an arbitrary network (e.g. a ST or a SN) and recursively construct the heights starting from the outlets with the aid of Eq. (2). From the derived landscape, a new network (in general different from the original one) can then be obtained by assuming that at each site the outward direction is along the steepest de- scent path and using Eq.(9). The noise term in Eq. (10) is mimicked by the unity term in Eq. (9). The proce- dure can then be iterated until self-consistency is finally achieved. The final configuration is, by definition, a sta- tionary state of Eq. (10). The convergence of the above procedure as a function

of the number of iterations is reported in Fig. 3 for various values of γ. Two remarks are in order. First we note that the dissipated energy as obtained from the above SC procedure does not have any physical meaning in the transient state. In other words the definition of “time” for the independent variable appearing in Fig. 3 should be considered only as a short-hand notation for “number of iterations”. Second it is apparent a the value γ = 0.5 seems to be the one with the slowest convergence ratio. This is probably due to the particular role played by the value γ = 0.5 as it will become clear shortly. Figs.4-8 depict typical patterns obtained on changing

the parameter γ in the interval [0, 1]. The effect of the value of γ on the aggregation pattern is evident. As γ increases from 0 to 1, single big rivers draining the entire basin in a snake-like form are less and less favored. The pattern for γ = 0 has a strong memory of the original initial tree. The overall effect of the SC procedure is to disfavor long meandering of the streams thus providing a self-affine character to the final tree. The main river be- comes rather straight for γ = 0.25 and the whole pattern appears to be more symmetric with respect to the case of γ = 0. A noteworthy feature is that this tendency is inverted as γ → 0.5 and returns back for γ = 0.75. As γ → 1 rivers become very directed as one would expect on the ground that this limit corresponds to the Schei- degger network [18]. The fact that γ = 0.5 most closely resembles real rivers is a reflection of the natural selec- tion by the erosional processes of this value of γ in terms of Eq. (2) as it has already been well documented in the literature (see e.g. [1]). In our numerical estimates of the exponents, for sim-

plicity, we shall restrict our attention to the half region [0, 0.5] . We also note that this is the region inaccessible to the analytical scheme of Ref. [17]. For a more accurate evaluation of the critical expo-

nents τ and ψ it proves convenient to introduce the inte- grated probabilities:

P (a, L) =

∫ a 0

da′p(a′, L) = a1−τF ( a

Lφ ) (11)

and

4

Π(l, L) =

∫ l 0

dl′pi(l′, L) = l1−ψG( l

Ldl ) (12)

where F (x) and G(x) are related to f(x) and g(x) de- fined in Eqs. (4) and (5), in an obvious way. An efficient way of computing the exponents is through the so-called “effective” (sometime also referred to as “running”) ex- ponents. In the present case they are defined as:

τ(a) = 1− ∂ logP (a, L)

∂ a (13)

and

ψ(l) = 1− ∂ logΠ(l, L)

∂ l . (14)

One then obtains an effective exponent for each value of the independent variable (a or l). Fig. 9 shows one typical result on a 256× 256 lattice.

We can divide the obtained values roughly in four regions. The first region (1 ≤ a < 10) corresponds to the region of no scaling. Small rivers belonging to the second region (10 ≤ a < 100) have an exponent close to the absolute minimum value τ = 1.5. This is consistent with the pic- ture of typical rivers (see Fig. 6) where small rivers dis- play a marked straightness similar to the one of the abso- lute minimum [18] and it means that they quickly assume configurations consistent with the absolute minimum (in- dependent of the initial conditions). Larger areas are as- sociated with larger rivers which have longer memory of the initial condition (ST in the present case). Hence the corresponding exponent is sensibly smaller (closer to the ST value 1.38). The last region corresponds to the finite size cut-off and must be discarded. After discarding the first and forth regions the ob-

tained values can then be grouped into local bins and a local average exponent can be associated with each of them. Statistical fluctuations within each box then yield an estimate of error bars. This provides our best esti- mate of the exponent for each value of L and a simple 1/L extrapolation is then carried out to extract the final values. This is depicted in Fig. 10 and the corresponding best estimates of this method are reported in Table I. An analogous procedure leads to the best estimates for ψ as reported again in Table I. Sizes and statistics are identical to those considered for the initial conditions and hence these simulations are rather time consuming. One can notice a weak dependence on γ for both τ and

ψ. One way of computing φ and dl is through the collapse plots of the probabilities P (a, L) and Π(l, L). However we have found that a satisfactory collapse can achieve only within a limited range of the appropriate variable (a/Lφ and l/Ldl). Hence we have opted for an alternative scheme hinging on the calculation of the following ratios:

M qa (L) ≡

⟨ aq+1

⟩ a

〈aq〉a ∼ Lφ q = 1, 2, ... (15)

where averages are over the probability densities p(a, L) (of the maximum river) and the L dependence is straight- forward [4]. A similar relationship holds for the lengths:

M ql (L) ≡

⟨ lq+1

⟩ l

〈lq〉l ∼ Ldl q = 1, 2, ... (16)

The results for q = 1, 2, 3 are reported in Table III for the exponent φ and in Table IV for the exponent dl. Sur- prisingly the values for φ are consistently larger that the space-filling value φ = 2 which is expected to be the up- per bound for this exponent. This is probably due to a deficiency of this procedure during a cross-over from a self-similar regime (the initial ST) to a self-affine pattern (the final configuration). Indeed we shall see later on that this feature is not present when one starts with a self-affine network (e.g. a SN) from the outset. A seem- ingly large value of dl is appearing probably due to the same reason. Overall, one can notice a rather weak (if any) dependence on γ and almost no dependence on q. Finally Hack’s exponent h was computed from its def-

inition Eq. (3) using an effective exponent method to be described below. Let us assume a power-law dependence for a generic

function as given by

f(x) ∼ xθ θ > 0 (17)

On integrating f(x) (between lower and upper limits, say x0 and x), we find

θ =

[ xf(x) − x0f(x0)

] ∫ x x0 dzf(z)

− 1. (18)

The effective exponent obtained using Eq.(18) in Eq.(3) with f(x) ≡ a, x ≡ l, and θ ≡ 1/h, can then be ana- lyzed with the same procedure (local average plus 1/L extrapolation) outlined before. We note that for a we have used an averaged value over all areas correspond- ing to the same length. The final results are reported in Table I and one can see that there is an overall good agreement with scaling laws.

C. OCN

The minimization of the energy functional Eq. (8) goes through an algorithm akin to the one exploited by Sun et. al [16]. It is based on the following steps:

1) An initial configuration (ST or SN) is generated and its dissipated energy is computed according to Eq. (8). By definition, this initial network has no loops.

2) A link is randomly selected and its local outflow is also randomly chosen.

5

3) This new configuration is tested for loop creation. If a loop has been created, the configuration is rejected and step 2) is repeated.

4) The energy of the new candidate configuration is com- puted. If it is smaller than the previous one, the new configuration is accepted, otherwise it is re- jected.

5) Steps 2)-4) are repeated until the energy does not change within a given tolerance.

The final configuration is regarded as a local minimum. This scheme is patterned after a standard Metropolis al- gorithm at zero temperature, and averages are over many different configurations (ranging from 500 at L = 32 to 100 at L = 256). Fig. 11 depicts the dissipated energy per unit of length as a function of the convergence “time” (i.e. the number of total iteration of the algorithm). The similarity with patterns obtained from the SC procedure is evident. Here, however, the typical convergence time is much longer than before. Critical exponents are computed with the same pre-

scription given in the previous subsection. In Fig. 12 we show the resulting 1/L extrapolation. The final best estimates are reported in Table I. Once more, a weak dependence on γ (τ increases as γ

increases) can be noticed. All other exponents are consis- tent with scaling relations Eq. (6). It is worth stressing that exponents are nearly consistent (at the edge of sta- tistical errors) with those previously obtained from the dynamical equation. Regarding the exponents φ and dl, they can be found in

Tables III and IV respectively. Here too the same feature, discussed in VB in connection with the exponent values of φ and dl, applies.

D. Independence of the initial condition

Our final task is a test of the sensitivity of critical ex- ponents to the initial conditions. To this aim we have changed the initial condition from a ST to a SN for both the dynamical equation and the OCN. The main differ- ence between the two initial conditions is that while the latter is a directed network, the former is not. Table II reports the comparison for the case γ = 0.5, and Fig. 13 depicts the network resulting from the SC scheme, for a typical final state. Despite the obvious memory of the initial network (see Fig. 2), the typical distance among big rivers is larger than the initial one. This in turn yields a higher value for H and hence non-trivial exponents. It is remarkable that although these patterns appear

to be considerably different from those obtained when starting with a ST, the two sets of critical exponents are nevertheless consistent within the statistical errors. As a final comment, we find in this case values of φ

and dl in very good agreement with scaling predictions. The results for different ratios q in the case γ = 0.5 for

both the SC and the OCN, along with the comparison with the corresponding values stemming from STs are reported in Tables V and VI respectively.

VI. CONCLUSIONS

In this paper we have addressed the issue of the exis- tence and robustness of the universality class associated with landscapes corresponding to local minimal energies. To this aim we have first extended previous studies

for both the SC solution of a Langevin equation and the OCNs variational methods. Higher sizes and statistics have been exploited in both cases and, (to our knowledge) for the first time, the most physical procedure of basing the statistics only on the river with largest flow, has been implemented. Second, we have monitored the dependence of critical

exponents on a parameter γ associated with the slope- area law in one case (SC) and with the dissipated energy in the other (OCN). Our results give compatible critical exponents between SC and OCN within the error bars, but a weak and similar dependence on γ appears in both models. Finally we have tested the stability of the obtained

results for both the SC and the OCN with respect to changes in the initial conditions. Although the obtained final patterns display a dependence on initial conditions, critical exponents appear to be insensitive to this depen- dence. As a by-product of our investigation, we have found

that the SC solution of the dynamical equation is a very powerful method to investigate river networks as it is ca- pable of providing useful informations on the stationary state in a simple and physical way. Another interesting point, from a numerical point of view, is that this pro- cedure typically achieves convergence much faster with respect to OCN scheme. In view of our results we can now summarize the argu-

ments favoring and disfavoring the appearance of a new universality class associated with local minima. As we mentioned in our discussion of Fig. 10 the typical evolu- tion of a network appears to depend on the considered length scale. Small rivers very quickly settle to their final state whereas much longer time is required to large rivers to forget their initial conditions. This is also reflected by the difficulty in collapsing the distribution probabilities of both a and l into a single plot for a reasonably ex- tended range of the corresponding variable. It is then possible that although only two universality classes (one associated with the initial condition and the other to the absolute minimum) are present, an intermediate univer- sality class sets on due to both the averaging over dif- ferent regimes and the difficulty of reaching the absolute minimum. On the other hand, in support of the existence of a new

universality class, we cite the fact that critical exponents

6

are different from those of both ST and SN – on start- ing with these initial conditions, critical exponents in the final configurations clearly deviate from their initial val- ues. Furthermore all exponents are found to be robust and obey to scaling relations summarized in Sec.III. Overall we believe that the evidence suggested by our

results favors this second possibility which has also been hinted into a different context [29]. In this respect the weak γ dependence of critical exponents remains unex- plained. It would be interesting to generalize the present cal-

culation in two aspects. For the dynamical equation, it would be instructive to tackle the problem of the ex- plicit solution of Eq. (1). This route has the advantage that Eq. (2) (the key relation between erosional process and network topology) can be then derived rather than assumed as in the self-consistent procedure. Similarly, a parallel calculation could also be implemented in the OCN framework, upon starting with the more general expression given in Eq. (7).

ACKNOWLEDGMENTS

It is a pleasure to thank Jayanth Banavar, Amos Maritan and Andrea Rinaldo for many enlighting dis- cussions, innumerable suggestions and a critical reading of the manuscript. Financial support from the Italian MURST (Ministero dell’Universita´ e della Ricerca Sci- entifica e Tecnologica) and INFM (Istituto Nazionale di Fisica della Materia) is gratefully acknowledged.

FIG. 1. A spanning tree on a 64× 64 square lattice with 8 nearest-neighbors, 4 in the Nord-South and East-West direc- tion and 4 along the diagonal directions. The outlet is on the West side (outflowing), on the East side there is an infinite wall, whereas there are periodic boundary conditions in the North-South direction. The thickness of the line at each point is proportional to the flow through that point. The seeming loops are just an artifact of the drawing.

FIG. 2. A Scheidegger network with L = 64. Boundary conditions are the same as in Fig.1. The directness of the network provides a privileged East-West direction.

7

0 500 1000 t

2500

3500

4500 E

/L

γ=0.00 γ=0.25 γ=0.50

FIG. 3. Dissipated energy per unit of length, as a function of the number of iterations (time t) of the self-consistent solu- tion of the dynamical equation, on starting with the spanning tree of Fig 1. The parameter γ is the value appearing in Eq. (2). The system size is L = 512 and each point of the curve is an average over all configurations which have gone at least that far in the number of iterations. All quantities reported in this and following figures are dimensionless.

FIG. 4. A typical network obtained as a final output of the self-consistent procedure described in Sec. VB. Here L = 64, γ = 0, and the initial condition is the spanning tree of Fig. 1.

FIG. 5. Same as in Fig. 4, with γ = 0.25.

FIG. 6. Same as in Fig. 4, with γ = 0.5.

8

FIG. 7. Same as in Fig. 4, with γ = 0.75.

FIG. 8. Same as in Fig. 4, with γ = 1.00.

0 200 400 600 800 1000 a

1.0

1.1

1.2

1.3

1.4

1.5

1.6

1.7

1.8

1.9

2.0

τ

FIG. 9. Effective exponent τ as a function of the area a in the stationary state of the dynamical equation with γ = 0.5. Here L = 256 and the initial condition is the spanning tree of Fig. 1.

0.00 0.01 0.02 0.03 0.04 1/L

1.0

1.1

1.2

1.3

1.4

1.5

1.6

1.7

1.8

1.9

2.0

τ

γ=0.00 γ=0.25 γ=0.50

FIG. 10. Finite size 1/L extrapolation for the value of τ obtained with the dynamical equation for the cases γ = 0, 0.25, 0.5. In all cases the initial condition is the span- ning tree of Fig. 1.

9

0 50000 100000 t

500

1000

1500

2000

2500 E

/L

γ=0.25 γ=0.50

FIG. 11. Dissipated energy per unit of length, as a function of the number of iterations in the OCN procedure. Here the size is L = 256 and again each point of the curves are averaged over all configurations which have reached that time t. The initial condition is the spanning tree of Fig. 1.

0.00 0.01 0.02 0.03 0.04 1/L

1.0

1.1

1.2

1.3

1.4

1.5

1.6

1.7

1.8

1.9

2.0

τ

γ=0.25 γ=0.50

FIG. 12. Finite size 1/L extrapolation for the value of τ obtained with the OCN (T = 0) for the cases γ = 0.25, 0.5. In both cases the initial condition is the spanning tree of Fig. 1.

FIG. 13. A typical 64 × 64 network obtained as a final output of the self-consistent procedure described in Sec. VB upon starting with a the Scheidegger network of Fig. 2.

TABLE I. Critical exponents τ , ψ and h as a function of γ for both the Self-Consistent (SC) procedure and the optimal channel networks (OCN) at zero temperature. The values in the parenthesis are those obtained from the computed value of τ and using the scaling relations. For γ = 0, the OCN exponents reported are the exact ones corresponding to the initial conditions (a spanning tree in the present case) since the energy is a constant.

γ τOCN τSC ψOCN ψSC hO 0.00 1.38 1.40± 0.01 1.6 1.69± 0.03 (1.67) 0.625 0.25 1.42± 0.01 1.42± 0.01 1.71± 0.07 (1.72) 1.68± 0.05 (1.72) 0.64± 0 0.50 1.44± 0.01 1.46± 0.01 1.8± 0.1 (1.79) 1.82± 0.05 (1.85) 0.61± 0

TABLE II. Critical exponents τ , ψ and h for γ = 0.5 as a function of the initial conditions. As in the text, ST and SN stand for Spanning Trees and Scheidegger Network respec- tively. SC and OCN have the same meaning as above. Values in parenthesis are scaling predictions.

τOCN τSC ψOCN ψSC hOCN ST 1.44 ± 0.01 1.46 ± 0.01 1.8 ± 0.1 (1.79) 1.82± 0.05 (1.85) 0.61± 0.03 SN 1.44 ± 0.03 1.43 ± 0.02 1.8 ± 0.2 (1.79) 1.7± 0.1 (1.75) 0.61± 0.03

10

TABLE III. Critical exponent φ stemming from both the SC procedure and the OCN scheme, as a function of the value of γ and of the value of the order q of ratios defined in Eqs. (13) and (14).

q γ = 0.00 (SC ) γ = 0.25 (SC ) γ = 0.50 (SC ) γ = 0.00 (OCN) γ = 0.25 (OCN) γ = 0.50 (OCN)

1 2.1 ± 0.1 2.2± 0.1 2.2 ± 0.1 - 2.0± 0.1 2.1± 0.1 2 2.1 ± 0.1 2.2± 0.1 2.2 ± 0.1 - 2.0± 0.1 2.1± 0.1 3 2.1 ± 0.1 2.2± 0.1 2.3 ± 0.1 - 2.0± 0.1 2.0± 0.1

TABLE IV. Critical exponent dl obtained from both the SC procedure and the OCN scheme, as a function of the value of γ and of the value of the order q of ratios (13) and (14) .

q γ = 0.00 (SC ) γ = 0.25 (SC ) γ = 0.50 (SC ) γ = 0.00 (OCN) γ = 0.25 (OCN) γ = 0.50 (OCN)

1 1.3 ± 0.1 1.3± 0.1 1.1 ± 0.1 - 1.2± 0.1 1.3± 0.1 2 1.3 ± 0.1 1.4± 0.1 1.2 ± 0.1 - 1.3± 0.1 1.3± 0.1 3 1.3 ± 0.1 1.4± 0.1 1.2 ± 0.1 - 1.3± 0.1 1.3± 0.1

TABLE V. Comparison between exponents φ as obtained from both the SC procedure and the OCN scheme, as a func- tion of the value of the initial conditions and of the value of the order q of ratios (13) and (14). As in the text, ST stands for Spanning Trees and SN for Scheidegger Network. Here γ = 0.5.

q φST (SC ) φSN (SC ) φST (OCN) φSN (OCN)

1 2.2 ± 0.1 1.8± 0.1 2.1± 0.1 1.7± 0.1 2 2.2 ± 0.1 1.8± 0.1 2.1± 0.1 1.7± 0.1 3 2.3 ± 0.1 1.8± 0.1 2.0± 0.1 1.7± 0.1

TABLE VI. Comparison between exponents dl resulting from both the SC procedure and the OCN scheme, as a func- tion of the value of the initial conditions and of the value of the order q of ratios (13) and (14). As in the text, ST stands for Spanning Trees and SN for Scheidegger Network. Here γ = 0.5.

q dlST (SC ) dlSN (SC ) dlST (OCN) dlSN (OCN)

1 1.1 ± 0.1 0.9± 0.1 1.3 ± 0.1 0.9± 0.1 2 1.2 ± 0.1 1.0± 0.1 1.3 ± 0.1 1.0± 0.1 3 1.2 ± 0.1 1.0± 0.1 1.3 ± 0.1 1.0± 0.1

[1] I. Rodrigues-Iturbe and A. Rinaldo, Fractal River Basins: Chance and Self-Organization (Cambridge Uni- versity Press, Great Britain, 1997).

[2] A. Rinaldo, I. Rodrigues-Iturbe, R. Rigon and R.L. Bras, Phys. Rev. Lett. 70, 822 (1993).

[3] P. Meakin, J. Feder, and T. Jossang, Physica A 176, 409 (1991).

[4] A. Maritan, A. Rinaldo, R. Rigon, A. Giacometti and I. Rodrigues-Iturbe, Phys. Rev. E 53, 1510 (1996).

[5] J.R. Banavar, F. Colaiori, A. Flammini, A. Giacometti, A. Maritan and A. Rinaldo, Phys. Rev. Lett. 78, 4522 (1997).

[6] K. Sinclair and R. C. Ball, Phys. Rev. Lett. 76, 3360 (1996).

[7] E. Somafai and L. M. Sander, Phys. Rev. E 56, R5 (1997).

[8] P. S. Dodds and D. H. Rothman, Phys. Rev. E 59, 4865 (1999).

[9] L.B.Leopold and W.B.Langbein, U.S. Geol. Surv. Prof. bf 500-A,1 (1962).

[10] M. Marsili, A. Maritan, F. Toigo and J. R. Banavar, Rev. Mod. Phys. 68, 963 (1996).

[11] R. Pastor-Satorras and D. H. Rothman, Phys. Rev. Lett. 19, 4349 (1998).

[12] F. Me´tiver, Phys. Rev. E 60, 5827 (1999) [13] A. Giacometti, A. Maritan, and J. R. Banavar, Phys.

Rev. Lett. 75, 577 (1995). [14] I. Rodriguez-Iturbe, A. Rinaldo, R. Rigon, R. L. Bras,

and E. Ijjasz-Vasquez, Water Resour. Res. 28, 1095 (1992); I. Rodriguez-Iturbe, A. Rinaldo, R. Rigon, R. L. Bras, and E. Ijjasz-Vasquez, Geophys. Res. Lett. 19, 889 (1992); A. Rinaldo et. al, Water Resour. Res. 28, 2183 (1992).

[15] S. Kirkpatrick, C.D. Gelatt, and M. P. Vecchi, Science 220, 671 (1983).

[16] T. Sun, P. Meakin and T. Jossang, Phys. Rev. E 49, 4865 (1994).

[17] A. Maritan, F. Colaiori, A. Flammini, M. Cieplak and J. R. Banavar, Science 272, 984 (1996).

[18] F. Colaiori, A. Flammini, A. Maritan and J. R. Banavar, Phys. Rev. E 55, 1298 (1997).

[19] The Peano basin is discussed e.g. in B.B. Mandelbrot The Fractal Geometry of Nature, (Freeman, New York 1983). Some topological properties in the river network framework have been addressed in A. Marani, R. Rigon and A. Rinaldo, Water Resour. Res. 27, 3041 (1991).

[20] We note that γ = 0 and γ = 1 correspond to known and exactly solvable models (spanning trees and Scheidegger network respectively) and that the only physically ac- ceptable values of γ lie between these two limiting cases.

[21] J. T. Hack, U.S. Geological Survey Professional Papers 294-B, 45 (1957).

[22] M. Cieplak, A. Giacometti, A. Maritan, A. Rinaldo, I. Rodrigues-Iturbe and J. R. Banavar, J. Stat. Phys. 91, 1 (1998).

11

[23] G. Caldarelli, A. Giacometti, A. Maritan, I.Rodrigues- Iturbe and A. Rinaldo, Phys. Rev. E 55, R4865 (1997).

[24] F. Y. Wu, Rev. Mod. Phys. 54, 235 (1982). [25] S. S. Manna, D. Dhar and S. N. Majumadar, Phys. Rev.

A 46, R4471 (1992). [26] S.S. Manna and B. Subramanian, Phys. Rev. Lett. 76,

3460 (1996). [27] M. Takayasu, H. Takayasu and G. Huber, J. Stat. Phys.

65, 725 (1991). [28] M. R. Swift, F. Colaiori, A. Flammini, A. Maritan, A.

Giacometti and J. R. Banavar, Phys. Rev. Lett. 79, 3278 (1997).

[29] B. Tadic, Phys. Rev. E 58, 168 (1998).

12

Comments