{'fraction_non_alphanumeric': 0.03894403268470262, 'fraction_numerical': 0.013854333080166566, 'mean_word_length': 4.141933746296795, 'pattern_counts': {'":': 0, '<': 5, '<?xml version=': 0, '>': 16, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 11, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
An Analytical and Numerical Study of Optimal Channel Networks

We analyze the Optimal Channel Network model for river networks using both analytical and numerical approaches. This is a lattice model in which a functional describing the dissipated energy is introduced and minimized in order to find the optimal configurations. The fractal character of river networks is reflected in the power law behaviour of various quantities characterising the morphology of the basin. In the context of a finite size scaling Ansatz, the exponents describing the power law behaviour are calculated exactly and show mean field behaviour, except for two limiting values of a parameter characterizing the dissipated energy, for which the system belongs to different universality classes. Two modified versions of the model, incorporating quenched disorder are considered: the first simulates heterogeneities in the local properties of the soil, the second considers the effects of a non-uniform rainfall. In the region of mean field behaviour, the model is shown to be robust to both kinds of perturbations. In the two limiting cases the random rainfall is still irrelevant, whereas the heterogeneity in the soil properties leads to new universality classes. Results of a numerical analysis of the model are reported that confirm and complement the theoretical analysis of the global minimum. The statistics of the local minima are found to more strongly resemble observational data on real rivers.


Introduction
Experimental observations on river networks have shown clear evidence of their fractal character. Data from many basins with different geological features have been analyzed, and have shown power law behaviour of the probability distributions of several quantities describing the morphology of the river basin [1,2,3,4] over a wide range of scales. Several statistical models have been proposed [5,6,7,8,9,10] but a complete theoretical understanding is still lacking. Recently a simple lattice model derived from an energy minimization principle has been proposed [11,12,13] that in spite of its simplicity, seems to reproduce many features of natural river networks.
Numerical investigations of the model have been performed [6,14,15,16] with different geometrical constraints on the form of the basin. Furthermore, the model has been analyzed within the framework of a finite size scaling Ansatz [17]. In the present paper the so called Optimal Channel Network model [11,12] is studied analytically and exact results are obtained. In addition, generalized models, taking into account the presence of quenched disorder are considered. Randomness is introduced in two different ways: one modelling the inhomogeneity of the soil and the other non-uniformity in the rainfall. Analytical results are given in these two cases, too. Results of numerical simulation in the last part of the paper, confirm and complement the analytic predictions.
In Section 2 we describe the lattice model and derive the scaling laws. The relationship between exponents are also derived. The exponents characterizing the power law distributions of drained areas and mainstream lengths are expressed in terms of a single independent exponent: the wandering exponent in the self-affine case and the fractal dimension for the self-similar basin. The section ends with the definition of the Optimal Channel Network model and with a short discussion of the underlying minimization principle. Section 3 is entirely devoted to a analytical study of the model on a simple fractal, the Peano basin. The solution is given exactly and is used to give bounds in the subsequent Section. The distributions of areas and lengths are evaluated and shown to exhibit the form predicted by the finite size scaling Ansatz. In Section 4, analytical results in the absence of disorder are derived. The model is shown to exhibit three distinct universality classes for different values of a parameter characterizing the dissipated energy. Heterogeneities in the soil properties and random rainfall are considered in the generalized models of Section 5. Analytical results for these cases are also deduced. Numerical results pertaining to the search for the global minimum of the dissipated energy with a simulated annealing algorithm are given in Section 6 and numerical results for the statistics of the local minima are given in section 7.

Definitions and Derivation of Scaling Laws
A river basin is described by a scalar field of elevations. Drainage directions are identified by steepest descent, i.e. by the largest local decrease of the elevation field. The presence of lakes has not been considered, i.e. from each point the water can flow to one of the nearest neighbours. This corresponds to have all lakes saturated. Within this context, a river network can be represented by an oriented spanning-tree over a two dimensional lattice of arbitrary size and shape, in which oriented links (one coming out from each site of the lattice) correspond to drainage directions. We will consider spanning trees rooted in a corner of a L × L square lattice (the root will correspond to the outlet). Site i is upstream respect to site j if there is an oriented path from i to j. To each site i of the lattice, we associate a local injection of mass r i (the average annual rainfall in the site i). The flow A i , also referred to as the accumulated area can thus be defined as the sum of the injections over all the points upstream of site i (site i included). The variables A i are thus related by where w i,j is 1 if site j is upstream with respect to site i and is a nearest neighbour of it and 0 otherwise. The local injection r i is commonly assumed to be homogeneous and identically equal to 1. In natural basins these areas can be investigated through an experimental analysis of digital elevation maps (DEM's) [2]. See figure 1a for an example.
The upstream length relative to a site is defined as the length of the stream obtained starting from the site and repeatedly moving in the upstream direction towards the nearest-neighbour with biggest area A (the one leading to the outlet is excluded, since it is a downstream site), until a source is reached i.e. a site with no incoming links (see Fig. 1b). If two or more equal areas are encountered, one is randomly selected. For a given tree, one may consider the probability distribution of the following quantities: for a lattice of given linear size L we will call p(a, L) the probability density of accumulated areas a and π(l, L) the probability density of the upstream lengths l. These represent the fraction of sites with area a or stream length l respectively. We will also consider the integrated probability distributions P (a, L), the probability to find an accumulated area bigger then a and Π(l, L), i.e. the probability to have a site with an upstream length bigger than l. Both these probability distributions, here defined in the simple case of the lattice model, were originally introduced to describe real rivers and experimentally found to scale as power laws leading to the formulation of a finite size scaling ansatz [17] p(a, L) = a −τ f a a C , where f (x) and g(x) are scaling functions accounting for finite size effects and a C and l C are the characteristic area and length respectively. The functions f (x) and g(x) are postulated to have the following properties: when x → ∞ they go to zero sufficiently fast to ensure normalization; when x → 0 they tend to a constant, to yield simple power law behaviour of the probability distributions in the large size limit. This also implies that τ and ψ are bigger than one. The characteristic area and length are postulated to scale as In river basins, anisotropies are always present due to a non-zero average slope of the landscape and the presence of gravity. Thus, one has to distinguish between a typical longitudinal length L and a typical perpendicular one L ⊥ (these two length are measured along the two principal axis of inertia), which scale as giving a C ∼ L 1+H , i.e. ϕ = 1 + H. H is known as the Hurst exponent and of course 0 ≤ H ≤ 1. In what follow, for the sake of simplicity, we consider basins of square shape; the above relations then, refers to dimension of a generic sub-basin inside the bigger one, whose dimensions are fixed from outside. The ϕ exponent thus corresponds to ϕ = 1 + H . The d l exponent, characterizing the typical length, can be assumed to be the fractal dimension of a stream (for fractal river networks, each rivulet going from any site to the outlet is a fractal with the same fractal dimension) , and is such that 1 ≤ d l ≤ 1 + H. The bounds correspond to a straight line and a space-filling behaviour. For self-affine river basins we expect d l = 1 and H < 1, whereas, when H = 1 then d l > 1 in the self-similar case.
For the same quantities, the integrated probability distributions can be analogously written: which follow from (2) and (3) with where sums over variable y have been replaced by integrals in large L limit. From the above definitions and the properties of f it easily follows that a n ∼ L (1+H)(n−τ +1) , for any n > τ −1, while a n ∼ const. if n < τ −1. Note that both a and l have a lower cutoff which is one. This equation, evaluated for n = 1 gives for the average area, The mean accumulated area a can be easily shown to be equal to the distance from a site to the outlet, averaged over all sites. In effect, in the sum over all the downstreams (the rivulets going from each site to the outlet), the number of times each bond (of unit length) appears exactly equals the accumulated area of the associated site. Thus summing over all A i is equivalent to a sum over all the downstream lengths: l downstream can be evaluated replacing in the sum the distance of each point from the outlet measured along the stream with the corresponding euclidean distance d(x) to the power d l : This fact is general and the argument we used does not need the knowledge of the downstream length distribution. This distribution however can be explicitly derived at least in the case of directed networks. We call directed those networks such that each oriented link has positive projection along the diagonal oriented towards the outlet. The reason to introduce this class of networks is that river basins often have a quasidirected character, due to the fact that they are typically grown on a slope which gives a preferred direction to the flow. Moreover directed trees are much more simple to handle analytically than "undirected" ones. For such trees consider the set of 2L diagonals orthogonal to the one passing through the outlet: the downstream length is the same for all the points on the same diagonal. Thus the number of points at a given distance to the outlet can be easily seen to be: number of point at distance l = l + 1 l = 1, ..L 2L + 1 − l l = L + 1, ..2L (15) thus the probability density for the downstream lengths has the form of a power law with exponent −1 times a scaling function of the argument l/L: The first moment of this distribution again gives eq. (14) with d l = 1 which is the expected result for the fractal dimension of a directed tree. This result, together with 14 suggest that downstream length distribution might have the scaling form: for the general case. Equations (13) and (14) lead to the following expression for the average area: From eq.(19) and eq.(12) we get the scaling relation A well known hydrological law, Hack's law [18], relates the length of the longest stream l in the drained area to the drained area a of the basin: The accepted value of h is h = 0.57 ± 0.06 [19,20,21], whose difference from the Euclidean value 0.5 lead to the first suggestion of the fractal nature of rivers [1]. From equations (4) and (5) it follows that Together with π and p one can define the conditional probabilityπ(l | a) of finding a main stream with length l in a basin with accumulated area a. The simplest scenario is that (21) still holds and [17]π(l | a) is a sharply peaked function of l with respect to a, i.e. there exists a well-defined constraint between lengths and areas, or more generally [22],π For the probability density π, p andπ the following consistency equation must hold: that gives, in the large L limit relating the exponents in the distribution of lengths and in the distribution of accumulated areas.
The scaling relations (20) and (26) can be expressed in a simpler form, observing that both τ and ψ depend on d l and H only in the combination d l (1+H) = h, where h is the parameter appearing in Hack's law (21). Thus The exponents characterizing the distributions of accumulated areas and upstream length are thus related by the simple expression: For self-affine river basins and all exponents can be expressed in terms of the Hurst exponent H, giving: while in the self-similar case yielding: Note that in both cases, τ ≤ 3/2. The equality holds only when H = d l = 1 which corresponds to the mean field situation [23]. Likewise, h ≥ 1/2.
A recently formulated lattice model [11,12,24] based on a minimization principle, seems to reproduce quite well the main characteristics of river networks. In this model, a rule for selecting particular configurations in the space of spanning trees is given. The "right" configurations, called Optimal Channel Networks (OCN) are obtained on minimizing a dissipated energy, written as: where Q i is the flow rate (the mean annual discharge) in the bond coming out from the site i, ∆z(i) is the height drop along the drainage direction, and k i characterizes the local soil properties such as the erodability. It will be taken to be equal to one for each bond for homogeneous river networks. Given a field of elevations, drainage directions are usually identified by steepest descent, i.e. by the largest downward gradient ∇z(i) of the scalar field z(i). This allows us to obtain another expression for the dissipated energy on adding the following considerations: i) in the case of uniform rainfall in time and space ii) experimental observations on rivers suggest an empirical relationship between the flow rate and the drop in elevation: with a numerical value around 0.5 for γ. Thus one obtains, apart from a multiplicative constant, the alternative expression of equation (36): which was first proposed by Rinaldo et al. [11,12,24], and will be analyzed further in this paper.

The Peano Basin
The Peano basin is a deterministic space filling fractal on which exact calculations can be carried out [25]. It has a spanning tree like structure not too dissimilar to that of real basins. The scaling laws for such a basin can be obtained exactly and the dissipated energy (37) can be estimated. The latter provides an upper bound for the minimum energy dissipated by an OCN and will be a crucial ingredient for the calculations to follow.
The Peano basin is obtained as follows: start with an oriented link; at step 2 such a link generates four new links, two resulting from the subdivision in half of the old link and preserving its orientation, and the other two having a common extreme in the middle point of the old link and both oriented toward it (see Fig. 2a). At each successive step, for each link four new oriented links are substituted in the way previously described. After T steps the fractal has N T = 4 T points (excluding the outlet) and it can be mapped on a square lattice of size L = 2 T with bonds connecting first and second neighbours to form a spanning tree. We can associate with each site i of the Peano basin (iterated until step T ) an area A i (T ) as in the previously defined lattice model of river networks. Let V T denote the set of distinct values assumed by the variable A i at step T . It can be easily checked that V T contains V T −1 and 2 T −1 new distinct values, appearing for the first time. Thus V T contains 2 T distinct numbers. Let us call A . = ∞ T =0 V T and a n the increasing sequence of numbers in A (the distinct values of A i that can be generated iterating the construction). For such a sequence, the following rule holds where the c k (n) are the coefficients of the binary expansion of n (n = k c k (n)2 k ).
In the construction described in Fig. 2b, denote by M T n the number of sites i with A i = a n at step T . The following recursion relation then holds: where t(a n ) is the first step in which an area with value a n appears, and is given by where [·] is the integer part. Solving (39) one gets and thus all the a n "born" at the same time step have the same probability p T (a n ) . = p(a n , L = 2 T ) = M T n /N T Then the integrated distribution of areas P (A i > a n , L = 2 T ) assumes a very simple expression for a n of the form 4 t (one can easily check from (38) that a 2 t −1 = 4 t ), and is given by having the form (7) with and and F (x) = 0 when x > 1.
Equation (42) is obtained on observing that P (A i > a = 4 t , L = 2 T ) = 2 T n=2 t p T (a n ) depends on n only through t(a n ), allowing one to replace the sum over n with a sum over the steps s. Moreover, for each step s > 0 there are 2 s−1 areas with the same t(a n ) = s. Thus Similarly, choosing l of the form l = 2 t and observing that at step T the sites with upstream length greater than or equal to 2 t are the ones in which the accumulated area exceeds 4 t , we find that which is of the form (8) with and In the following Section we will need some estimates of the mean value of A γ i for a Peano basin of size L = 2 t .
p T (a n )a n γ .
From the expression for a n , it follows that (t(a n ) is the "birth-time" for a n ), giving where Performing the summation one gets in the large size limit From equation (51), one gets which will be essential to obtain an energy bound for the lattice OCN model. The scaling exponents for the Peano basin can also be obtained by a renormalization group argument. Let us consider, for example, the scaling of the accumulated areas. The self-similar structure of the Peano basin suggests a natural "decimation" procedure. The idea is the following: consider the equations relating areas at time step T ; then, eliminate the variables related to the sites that are not present at time step T − 1. This leads to an effective equation describing the same physics on a tree scaled down by a factor 2. For the sake of simplicity let us consider the Peano basin at the second step of iteration. In Fig. 3 let A (2) n denote the variables related to sites that are present at step T = 1 and B (2) n denote the ones that will be eliminated by decimation. The upper label refers to the step. In what follows it will be useful to write the equations in terms ofÃ (T ) n = A (T ) n − 1 andB (T ) n = B (T ) n − 1. The areas at step T = 2 are related by:Ã 1 =Ã Elimination of theB (2) n leads toÃ At time step T = 1 the relation between areas is straightforward: Equations (56) and (57) are the same if i.e.
(A (T +1) Denoting by n (T +1) (a) the number of sites with area greater then a at step T + 1, one can easily observe that the number of decimated sites with A > a is half of the total number of sites with A > a n (T +1) (a) = 2 · n (T ) (a/4) .
Noting that the total number of sites at step T , N T = 4 T , the integrated probability P (A (T +1) n > a) = n (T ) (a) N T is given by The general solution of equation (61) is of the form apart from a superimposed oscillating term given by a periodic function of log(a) with period 1/4 [25]. From equations (61,62,63) The same argument can be repeated for the distribution of stream lengths, recovering the ψ exponent.

Analytical Results: homogeneous case
We now proceed to an analysis of the characteristics of the global minimum of the functional E for γ in the range [0, 1] for the homogeneous case. Let us consider first the two limiting cases γ = 0 and γ = 1. If we call l i the weighted length of the stream connecting the i-th site to the outlet (calculated assigning to each bond a length k i ), it is straightforward to show that In effect, denoting with DS(i) (US(i)) the set of points downstream (upstream) with respect to the point i and observing that A(i) equals the number of points in the set US(i) one gets: The minimization of the energy dissipation for γ = 1 thus corresponds to the minimization of the weighted path connecting every site to the outlet independent of each other.
The γ = 0 case, on the other hand, corresponds to the minimization of the total weighted length of the tree: In the homogeneous case k i = 1, ∀i, which leads to a high degeneracy for both γ = 0 and 1. Indeed for γ = 0, each configuration has the same energy (each spanning tree on a L × L square lattice has L 2 − 1 links). For γ = 1, the minimum of the energy is realized on a large subclass of spanning trees, namely all the directed trees, in which each link has a positive projection along the diagonal oriented towards the outlet.
For the values of γ ∈ (0, 1) there is a competition between both mechanisms breaking the degeneracy and making the search for the global minimum a less trivial problem. The γ = 1 case gives a minimum energy E ∼ L 3 . This can be derived observing that all points on a diagonal orthogonal to the one passing through the outlet have the same distance from the outlet. Then: E = L−1 k=1 k(k + 1) + 2L−2 k=L k(2L − 1 − k) = L 2 (L − 1) ∼ L 3 . Thus the value of the energy functional is the same for each directed network and corresponds to the Scheiddeger model of river networks [5] -all directed trees are equally probable. Such a model can be mapped into a model of mass aggregation with injection that has been exactly solved by Takayasu et al. [26,27]. The corresponding exponents are: These exponents follow easily from the result H = 1/2 and from our scaling solution of sec. 2. The result H = 1/2 can be deduced with a simple argument: since all directed trees are equally probable, each stream behaves like a single random walk in the direction perpendicular to the diagonal through the outlet. This implies that its perpendicular wandering is L ⊥ ∼ L 1/2 . Comparing with equation (6) one gets H = 1/2 . The γ = 0 case leads to the same energy E ∼ L 2 for each network, thus reducing to the problem of random two dimensional spanning trees, whose geometrical properties have been calculated on a square lattice [28,29]. The results, in our notation, are: Both (67) and (68) are consistent with the scaling relations (27), (28) and (29).
We now extend our analysis to the whole range γ ∈ [0, 1]. We will rigorously show [30] that, in the thermodynamic limit, the global minimum in the space S of all the spanning trees of the functional E(γ, T ) = i A i (T ) γ scales as for all γ ∈ [0, 1] .
Since E(γ, T ) is an increasing function of γ and it is equal to L 2 for γ = 0, for γ ≥ 0 it is obvious that We now observe that the sum over all the sites can be performed in two steps: where D n are the lines orthogonal to the diagonal passing through the outlet, that we will enumerate starting from the corner farthest from the outlet (see Fig. 4). For directed spanning trees one can observe that the sum of the areas in a given line D n is independent of the particular tree, and is: where S d (0) . = 0. This quantity only increases on generalizing to generic "undirected" trees and thus for an arbitrary spanning tree, Let us observe that for k = 0, ..., (L − 1) : making it convenient to perform the summation in (71) over "pairs" of lines. To get a lower bound for E we need a further inequality: for every set Γ that follows easily from Schwartz inequality, since A i ≥ 1 and 0 ≤ γ ≤ 1. Now, using (71), (74) and (75) we can write whereD n = D (2L−1−n) . The equality in the last inequality holds for directed networks. We can thus write Equation (77) together with the (70) gives the lower bound that holds for every T ∈ S and thus also for the minimum over T . Using the results of the previous Section, we can exhibit a tree on which the bound is realized. In effect, the Peano network can be mapped on a square lattice only considering links between first and second nearest neighbours, but (78) can be analogously obtained for such a lattice on rearranging in an opportune way the summation in equation (71). If we call T P the spanning tree given by the Peano basin we know from equation (54) that, except for logarithmic corrections for γ = 1/2, E(γ, T P ) ∼ max(L 2 , L 1+2γ ). Thus and we get equation (69). We now proceed to the calculation of the scaling exponents. For a directed path, from equation (65), For generical "undirected" networks, let us write as in (19) a ∼ L d l (81) (d l could be be somewhat bigger than one if one assumes a "quasi-directed" behaviour). Using eq. (11) with n = γ and the above result on the scaling of energy: E = L 2 a γ ∼ L 1+2γ , one gets holding for γ > τ −1 . Equations (82) together with the scaling relation (20) can be solved with respect to τ and H, and gives, for γ > 1/2 (the constraint γ > τ −1 become γ > 1/2, independently of d l ). Thus, if H ≤ 1, for any γ < 1 d l = 1, yielding for γ ∈ (1/2, 1). The exponents are the same as in the mean field theory [26,27]  When 0 < γ < 1/2, all we can say is that τ > 1 + γ. However if d l = 1 for any γ ∈ (0, 1) (i.e. the Optimal Channel is quasi − directed) then one might expect that H = 1 for all values of γ and thus eq. (84) holds in the whole range γ ∈ (0, 1). This prediction is confirmed by results of numerical simulations.

Analytical Results: heterogeneous case
In this Section we focus our attention on the case in which some sort of quenched disorder is present in the basin. Two cases have been analyzed: i) random bonds, modelling heterogeneity in the local properties of the soil; ii) random injection, modelling non-uniformity in the rainfall.
In the first case, we will show that the energy can be bounded from above with the corresponding one in the absence of disorder. This will give in the large L limit an upper bound for the τ exponent that we will see is realized for all the γ ∈ (1/2, 1). In the case of random rainfall, we will show that this kind of randomness does not affect the scaling behavior of the dissipated energy in the large L limit. All the analytical results found in Section 4 for the homogeneous case, being based on the energy estimate in the thermodynamic limit, can thus be extended to that case, giving the same values for the exponents.
In the case of random bonds, we associate with each bond of the L × L square lattice (2L(L − 1) bonds) a quenched random variable k b , arbitrarily distributed such that k b = 1. The label b ranges over all the bonds of the lattice. The 2L(L − 1) variables are chosen independent of each other and identically distributed.
In what follows b (i, T ) will denote the label associated with the bond coming out from the site i on the tree T . Let T * (γ) denote one of the trees on which the minimum of the energy E(γ, T ) is realized in the homogeneous case for a given value of γ, and S the set of all the spanning trees: Denoting by E D (γ) the minimum of the energy in the heterogeneous case, averaged over the quenched disorder, the following inequality holds: The energy in the presence of this kind of disorder is thus bounded from above by the energy in the absence of disorder for any value of the γ parameter. In the large size limit this result gives bounds on the scaling exponents. Equation (11) evaluated for n = γ gives which holds for any γ > τ D −1 = 1−h D , thus at least for any γ > 1/2 since h D ≥ 1/2.
Here and in what follows variables with the D index refer to the random bond case. Equation (87), compared with equation (69) leads to In the case of self-affine behaviour, using eq. (20), equation (88) gives: As before, for the γ = 1 case the above inequalities are useless. In the case of self similar behaviour, inequalities for the fractal dimension d l,D can analogously be deduced: Because, in general 0 ≤ H ≤ 1 and 1 ≤ d l ≤ 1 + H , equations (89) and (90) together yield for γ ∈ (1/2, 1): H D = 1, d l,D = 1, i.e. the mean field values. The γ = 1 case has been exactly solved [23], and it can be mapped into a problem of a directed polymer in a random medium or equivalently a domain wall in a disordered two dimensional ferromagnet [32,33,34], for which an exact solution is known.
The corresponding values for the exponents are: Disorder can be introduced in the system in another way, replacing the constant injection in each site of the lattice with a random, quenched, local injection, i.e. a spatial inhomogeneity in the rainfall r i in eq. (1). In order to do that, one can associate with each site i of the lattice a random variable r i . The variables are chosen to be independent of each other, identically distributed and with mean r i = 1 . The accumulated areas must then satisfy, as in (1) in such a way that The minimum of the energy averaged over the "random-rainfall" will be denoted by E rr (γ), and for a given value of γ is given by where S denotes the set of all spanning trees T and {r j } denotes the whole set of random variables. As in (85) T * will be one of the trees for which the minimum of the energy is realized in the absence of randomness in the rainfall and for a given value of γ. Then Thus In this case, it is possible to bound the energy also from below with an argument analogous to that used in Section 4 for the homogeneous case. The detailed calculation is given in the appendix. Thus one can conclude that and all the results of Section 4 hold.

Numerical Results: global minimum
We have carried out extensive numerical investigations of OCN along two avenues: i) the search for the global minimum with a Metropolis algorithm for γ = 1/2.
ii) the statistics over local minima for γ = 1/2; strikingly these yield consistent but different values for the scaling exponents. By a local minimum, we mean a configuration (a spanning tree) of the network such that no link can be changed without increasing the energy. The global minimum is of course also a local minimum; but in the two cases we found different statistics, that is suggestive of a very rich structure of the energy-landscape. We will postpone the results concerning the latter subject to the next Section, focusing our attention only on the scaling properties of the global minimum.
In the computer simulations we considered a square lattice with all sites on one side allowed to be an outlet fot the network. Periodic boundary conditions were chosen in the other direction. Once the simulation has been performed over the whole lattice, the basin with the biggest drained area is selected and only sites contained therein are used to calculate statistical quantities. Multiple outlets are allowed in order to minimize finite-size effects. The rainfall is assumed to be uniform over the whole lattice. The optimization method used is simulated annealing, in which a parameter T analogous to the temperature is introduced and lowered during the simulation. For each T value the system is relaxed in the following way: a new allowed configuration "near" the initial one ("near" means one differing from the previous one only in one link) is randomly chosen; the dissipation energy of the new configuration is calculated and compared with the value of the old one. The new configuration is accepted with probability 1 if ∆E is negative, and with probability exp[− ∆E T ] otherwise. In short a sketch of the algorithm is: i) Generation of a random initial configuration: Starting from a given, fixed tree we obtain a random initial configuration running steps ii) and iii) several times. In this way only changes preserving the spanning loopless structure are accepted.
ii) Random changes of the configuration: We select randomly one site i and then again randomly one of the "free" (not yet belonging to the tree) links connected with i, if any.
iii) Geometrical controls: The absence of a loop in the new configuration is checked. If a loop is present, step ii) is repeated. In each cycle, the T parameter is decreased with the following rule: at the n-th cycle T is given by T (n) = α n T (0), where α is a parameter very close to 1 (we choose α = 0.986) and T (0) is a suitable chosen constant.
After step i), steps ii)-v) are repeated many times, say N. Then we go to step vi) in which the T parameter is lowered and the entire algorithm (except step i) ) is repeated. The annealing process stops when T reaches very low values (≈ 10 −4 ).
The simulations have been repeated varying the initial configuration for a size L = 128. The statistical quantities do not depend on the initial data. The integrated probability distributions for the accumulated areas and mainstream lengths averaged over 10 trials are shown in Fig.5 and Fig.6 and give τ = 1.50±0.02 and ψ = 2.00±0.02 where the error is estimated as the root mean square root over the ten trials. These results are in perfect agreement with equation (29) and confirm that the analytical results hold when γ = 1/2.

Numerical Results: local minima
Homogeneous OCN yield results in good agreement with experimental data on rivers when statistics based on local minima is calculated. This suggested [35] the concept of feasible optimality according to which nature is "unable" to reach the true ground state when complex systems are involved in optimization problems. The optimization just stops when one of the local minima is reached. Within this framework, the scaling properties of real rivers should be reproduced considering the ensemble of local minima.
We have performed the search for local minima with an algorithm equivalent to a T = 0 Metropolis scheme, i.e. one in which new configurations are accepted only if energetically favourable. The simulation has been repeated 40 times, starting with different, randomly chosen, initial data and varying the size L of the system: L = 32, 64, 128 .
The values obtained for the characteristic exponents of the probability distribution considered above, τ and ψ are in very good agreement with experimental data. The distributions obtained starting from different initial condition do not substantially differ from one another, all yielding the same value for the exponents.
The statistical quantities thus seem to be robust with respect to variations of initial conditions, showing the self-averaging of the scaling parameters.
Results obtained averaging over 40 local minima are shown in Figs.7 and 8 and give τ = 1.45 ± 0.02 and ψ = 1.82 ± 0.02. The results are consistent with the scaling relations. A collapse test has been done to verify the consistency of numerical values of the exponents with the finite size scaling hypothesis (Fig.9).

Summary and Conclusions
In this paper, we have studied the Optimal Channel Network for the drainage basin of a river. Within the framework of the finite size scaling hypothesis for the distributions of accumulated areas and mainstream lengths, we deduced the exact scaling behaviour for the tree (drainage network) for which the absolute minimum of dissipated energy is realized in both the homogeneous case and in the presence of randomness. The scaling exponents in the homogeneous case are found to be the mean field ones and differ from those measured in real rivers. Numerical results were obtained both for the statistics of the global minimum (confirming the analytical predictions) and for local minima. The statistical exponents characterizing the local minima definitely differ from the mean field ones. They seems to be in a new distinct universality class and in agreement with data from real rivers. This suggests that real rivers, during their evolution, do not visit all of configuration space but are soon trapped in a metastable state, i.e. a local minima of the dissipated energy. The authors wish to thank A. Rinaldo and M. Cieplak for many enlightening discussions. This work was supported by NASA, NATO, NSF, The Petroleum Research Fund administered by the American Chemical Society and The Center for Academic Computing at Penn State.

Appendix
Using the notation of Section 4 we denote by D n the lines orthogonal to the diagonal passing through the outlet, enumerated starting from the corner farthest from the outlet with 0. It will be useful to associate each line D n with n ≤ L − 1 with aD n : In what follows we will choose the {r i } to be independent random variables with values in the interval [0, 2] and such that r i = 1. We then proceed to the first step of the proof as in Section 4, observing that the sum over all the sites in equation (94) defining the energy can be performed in two steps: where in the last step the equality holds only for directed trees. The last expression can be written in a more convenient form introducing Then where the last inequality follows on observing that Thus, for r i ≤ 2, the argument between the brackets is less than or equal to one. Furthermore, for any x and γ ∈ [0, 1], x γ ≥ x. Equation (101) together with equation (96) of Section 5 gives and thus [31] OCN networks on a square lattice have also been studied by Sun, Meakin and Jøssang [14]. They allowed each site of the border to be an outlet and considered the distribution of accumulated areas of a site belonging to the border (and not for a generic site of the basin as we have been doing in the present paper). They found a power law behaviour for such a distribution: P ∂ (a) ∼ a −τ ∂ with τ ∂ numerically estimated to be ≃ 1.51. If finite size scaling is invoked for such a distribution [Meakin, Feder and Jøssang [6]], P ∂ (a, L) = a −τ ∂ f ∂ a L 1+H then, as in (12), one finds: As noted in [6] the mean accumulated area for a border site in a basin is just the number of sites in the bulk divided by the number of sites in the border. Denoting by D the fractal dimension of the border, a ∂ L D ∼ L (1+H) , τ ∂ = 1 + D 1+H . In the case of compact basins (H = 1), τ ∂ = 1 + D 2 , ( τ ∂ = 3/2 for a square ). Note that τ ∂ is independent of the class of spanning tree chosen to drain the basin and it is a priori different from the τ we introduced to characterize the distribution of accumulated areas in the bulk, which we found to be for H = 1, τ = 2 − d l 2 .