Exact solution for the self-organized critical rainfall model

Exact solution for the self-organized critical rainfall model

R. F. S. Andrade

Instituto de Física, Universidade Federal da Bahia Campus da Federação, 40210-340, Salvador, Bahia, Brazil

ABSTRACT

This work presents an analytical investigation for a Self-Organized Criticality abelian model that describes basic properties of rainfall phenomena. The knowledge of the exact solution for the probability that a site topples when mass is added to any other site of the lattice leads to a large number properties of the model, including the exponent of the power law that describes presence of the events as function of their magnitude. It is shown that the model belongs to the same universality class of a first model proposed by Dhar and Ramaswamy (DR). However, for finite size lattices, it is found that its exponent is larger than that one for the DR model.

I Introduction

Abelian sand pile models (ASM) [1] became quite important for the understanding of basic properties of Self-Organized Criticality (SOC) [2, 3]. They satisfy the remarkable property that the final state of the model, after subsequent addition of grains to any two sites, is independent of the order in which the grains have been added.

The first analyses by Dhar and Ramaswamy [4] considered a model (DR) based on a critical height criterion for toppling, i.e., a site becomes unstable and topples when its amount of mass exceeds a threshold value m_th. In this model, a site topples to exactly d distinct neighbors. This is a directed and deterministic model as, i) the toppling rules breaks the isotropy of the lattice and avalanches develop along one direction; and ii) a deterministic rule indicates the fixed number of grains that each site receives when an unstable site topples. These two properties impose the condition that each site topples only once during an avalanche.

Abelian models constitute a special set where analytical investigation has lead to exact results, what are still rare subjects in the SOC landscape. More recently, a number of works have focused on certain variants of abelian models: they include the random distribution of toppling grains onto a restricted number of neighboring sites [6-8], and the complete toppling of all grains from an unstable site [9-11]. The first modification does keep the models in the abelian class, and an exact solution for the probability distribution function of events (PDF) has been discussed. The properties of such models are distinct from those with deterministic toppling rules. On the other hand, the second modification breaks the commutativity of the models, that seem to be non-integrable.

Other SOC models, that were initially proposed based on the so-called gradient condition for toppling, as the Bak-Tang-Wiesenfeld (BTW) model [2], also belong to this class, provided the variables are conveniently re-interpreted. But, as far as they allow for multiple toppling for a single avalanche, they are not exactly integrable.

The exact solution for the DR model is based on the evaluation of the probability distribution function P(s₁;s₀). It measures the probability that a site s₁ topples when one grain is added in a site s₀, and takes into account all configurations of the lattice which satisfies this requirement. The knowledge of P(s₁;s₀) opens the door to a large number of properties of the model, including the function r(f), that measures the relative number of avalanches r of a given size as function of the avalanche size f.

The ASM class includes a much larger number of models, some of which can be exactly integrable. A model suggested by a very crude description of drop avalanche inside a 2-dimensional cloud, the abelian rainfall model [5], is quite close to the original DR model. However, it allows for the presence of holes inside an avalanche cluster, what is found in the DR models only for dimension d > 3. Although the numerical simulations suggest that this model belongs to the same universality class as the d = 2 DR model, its analytical solution is still missing. Such solution can clear out whether the presence of holes inside an avalanche cluster changes or not the critical behavior of the model. The purpose of this work is to discuss a full analytical solution for P(s₁;s₀) of this model, and to prove that, in the infinite lattice limit, r(f) is indeed described by the same exponent as in the case of the DR model. We also analyze in which extent, results for finite size lattices for the two models differ from each other.

The rest of the work is so organized: in the Section II we present the model, writing the equations for P in terms of two sub-lattices; the corresponding solutions are obtained in the Section III. Section IV discusses the relation between P and r, and this function is evaluated with the help of the three sites probability function P₃. Section V brings results for finite size lattices, as well as the analytical solution in the infinite lattice limit. Section VI closes the work with concluding remarks.

II The two sub-lattices model formulation

The model for rainfall is defined on a N × M square lattice, whose axes are oriented along horizontal and vertical (downward) directions [5]. Periodic boundary conditions are imposed on the horizontal direction. Each site s_i on the lattice, labeled by the indices i = (j,k), is intended to describe a condensation nucleus, around which vapor condenses, leading to cloud droplets with liquid water content measured by the variable m_i. When a water ''grain'' is randomly added to a site i, representing the growth of droplet by mass aggregation due vapor condensation, it topples if the value for m_i > m_th. This describes the downward motion when the gravitation surmounts the drag force. It is followed by collision, coalescence and drop break-up. The mass of site i and its neighbors in the lower row are updated according to the rule:

where w + u + e = 1.

Let us consider a site s_i and its corresponding P(s_i;s₀) = P(j,k;s₀). Then, as the site s_i receives grains from its neighbors in the upward direction, we can relate P(s_i;s₀) with the same functions evaluated at these sites, i.e.,

Now we split the original lattice into two square sub-lattices, A and B, as indicated in Fig. 1. The axes of both sub-lattices are obtained from the old ones by a rotation of p/4 followed by an inversion of the new X axis. In this reference frame, the sites are re-labeled according to:

We observe that the row index k is expressed in terms of the new labels by + .

The equation (2.3) is rewritten in terms of two probability functions P_A and P_B, which are the restrictions of P to each of the two sub-lattices:

For the sake of simpler notation, we have let [ , ] ® [j, k] and j + k = t in (2.4) and in all expressions from now on.

If we set u = 0 in (2.4), the equations for the sub-lattices decouple and become equal to that one used to describe the DR model in two dimensions. This confirms that the two models are very close, but the solutions of the equation for general u may lead to other type of behavior.

III Solution of the equation for P

A particular solution to (2.4) can be evaluated after a series of steps. So let us define

and

In terms of these new functions, the equation (2.4) becomes

Bearing in mind the typical solutions for the DR model, we look for solutions in which the functions Z_j,k become separable, i.e.,

and use the following ansatz

Inserting (3.4) and (3.5) into (3.3), and imposing that the coefficients of the different powers of j and k must be equal on the different sides of the resulting expressions, we obtain several relations between A,B,b, z, h, x and the parameters e,w and u. The evaluation of the complete set of roots to this equation is a difficult task, as the unknowns are present in several factorials. However, the following particular solution, valid when e = w = u/2 = 1/m_th, has been found by inspection:

It leads to the explicit form

As expected, this solution satisfies

where t₀ = j₀ + k₀.

The probabilities P(s₁;s₀) and P(s₂;s₀), for s₁¹ s₂, are not independent, as there are many configurations that lead to toppling from both s₁ and s₂. In the sum of P(s;s₀) over all sites with t constant, the configurations that lead to q toppling sites appear exactly q times. The sums in (8) are equivalent to summing over individual configurations multiplied by the number of toppling sites. Thus, to the average number of sites that topple when one grain is added at s₀. As m_th = 4, this indicates that the average of one grain falls from any row to each added grain, and reflects the mass conservation property of the model.

It is worth pointing out that a general analytical solution, for any values for e,u and w, may not exist in terms of the ansatz (3.5), the form of which is too restrictive. However, solutions for other values of these parameters, though not expressed in a closed compact form, can indeed be found [12]. The present solution must capture the essential global properties of this model, as different choices for the parameters do not change its symmetry properties. This can also be seen in the results of many different numerical simulations.

IV The distribution of events r (f)

As the function P reflects the constant flux property, it does not directly measure the size f of avalanches, that is defined by the number of sites that topple. To obtain the function r(f) we first relate f to the number of rows that an avalanche reaches. To this purpose, we have to evaluate N(t;t₀), the average number of events in which the row t is active, i.e., in which at least one site in the row t topples when a grain is added at row t₀. It is easy to observe that N(t;t₀) = N(t – t₀;0), so that we simplify the notation to N(t), assuming that the grain is added at t₀ = 0. Thus, the probability (t) that one avalanche exceeds the row t is N(t)/C(t), where C(t) counts the number of distinct configurations of the lattice until row t, i.e., C(t) = , with S(t) = (t + 1)² indicating the total number of sites that can be reached by an avalanche. If (t), a decreasing function of t, is characterized by some exponent a, the average mass m(t) that topples from t must increase with t ^a , in order to keep the constant average flux. Then we obtain for the integrated mass f(t) ~ t ^a+1, so that after inverting this relation to t ~ f , we can express (f) º (t(f)) ~ f ^– . As (f) expresses the probability that an avalanche exceeds the total number of sites f, it is related to r(f) by

The behavior of the function r(f) depends essentially on the behavior of N(t). It can be expressed by

where n_t() counts the number of configurations in which sites topple in the row t. The behavior of n_t() depends on the values of the same function for t – 1, so that a general equation can be written as

Due to the presence of holes in the avalanche clusters, the form of the coefficients c_t(,^¢) for the rainfall model becomes very cumbersome as the value of t increases. However, for the DR model in two dimensions (u = 0 in (2.4) and m_th = 2), c_t(,>^¢) has the very simple form:

Such simple form for c_t(,^¢) is related to the fact that the avalanches for the DR model in two dimensions have no holes. Equation (4.3) with the coefficients given by (4.4) can be numerically integrated, so that the behavior for (t) and the value for a is easily obtained.

A possible way to side-step the difficult to evaluate (4.3) for the present model, is to consider the square average flux over the constant row t. It is obtained by summing, over all lattice configurations until that row, the square of the number of sites that topple from the row t for each particular configuration. It can also be regarded as the average flux taken with the help of other distribution that favors those configurations with larger number of toppling. In comparison to the constant flux (eq. (3.8)), this quantity over-weights the configurations with larger number of toppling sites, so that it increases with t. Its dependence on t follows the same rule as for m(t), so that it leads to the exponent a.

To evaluate the square flux over a given row t, we will consider the three point probability function P₃(s₁, s₂; s₀), where both s₁ and s₂ belong to the that row. Of course, we must have P(s₁; s₀) = P₃(s₁, s₁; s₀). This function essentially counts the number distinct of configurations until that row, where both s₁ and s₂ topple. Following the same arguments to identify (3.8) with the average flux, it is possible to see that, in the sum P₃(s₁, s₂; s₀), a given configuration with q toppling sites appear exactly q² times. As a consequence, this sum equals the square average flux, and should increase with t ^a .

The equation for P₃(s₁, s₂; s₀), that looks much the same as (2.4), now depends on two indices Z and Z^¢. For instance, taking Z = Z^¢ = A, we obtain:

where we dropped the explicit indication of the site s₀ for a simpler notation. The equations corresponding to the three other possible choices of Z, Z^¢ are very similar to (4.5).

The solution to the equations for P₃ can be obtained from the following ansatz:

where y – s₀ must be read as the difference between two vectors, and the sum over y spans all sites on the rows t^¢ such that t₀< t^¢ < t . Inserting (4.6) into (4.5) and the three other equivalent equations, and demanding that P(s₁; s₀) = P₃(s₁, s₁; s₀), leads to a set of equations that allows, recurrently, for the evaluation of F(y) at all different sites.

For the current purpose, however, it is sufficient to evaluate

where the prime indicates that the sums are restricted to the sites s₁ and s₂ over a given row t. Then we observe that

so that (4.7) becomes:

Noting that the function P_Z(s₁; y) depends only on the relative position of s₁ and y, we can re-label all positions in the above equation, setting t₀ = 0, so that (4.9) becomes

Defining (t) = F(y), and using (3.8) we obtain

To evaluate (t) let us consider a restricted version of (4.7), where the sums are taken for s₁ = s₂. In this case, (4.10) becomes

We define

so that (4.12) reduces to

This is a very compact way of writing the equations for F(t). It allows for a formal solution F(t) = F(t – 1) – F(0)d(t) that can be easily computed, with

V Results

The square flux has been evaluated for several values of t with the help of (4.11), as shown in a logarithm plot in the Fig. 2. For the purpose of comparison, we also draw, in this figure, results for the DR model: the square flux, that has been evaluated by a similar procedure as described above, and the function m(t), that has been evaluated with the help of eqs. (4.2-4.4). The curves for the DR model are parallel in the limit of large values of t, what confirms that both the square flux and m(t) depends asymptotically on t with the same exponent a. However, it is possible to note that the curve for the rainfall model is somewhat steeper.

For a more accurate analysis we show, in the Fig. 3, local values for the exponent (t) obtained from the square flux data for both DR and rainfall models. The evaluation of the corresponding a(t) proceeds via a least square evaluation on a local window around t. The curves were drawn for a window encompassing 10 points, but the results depends only very weakly on the window width. We see that the exponent for the DR model converges to the asymptotic value 4/3 obtained by Dhar from the lower side, according to power law |4/3 – t| ^–1. On the other hand, the exponent for the rainfall model is always greater than 4/3 and decays to this value with a lower exponent ~ -0.47.

As the crucial difference between the two models lye in the presence of holes inside an avalanche cluster, the results show that they cause the rainfall model to have a smaller probability of larger events in comparison to the DR model, specially for small lattices. In the thermodynamic limit, however, the two models belong to the same universality class.

A confirmation of the above result can be obtained by exploring the properties of the analytical solution for the square flux. We see from (4.13) that the contributions to Q(t) exhibit a sharp maximum at = t. Thus, we estimate

where D(t) indicates the width of the peak around = t. Making use of the Stirling approximation, we can easily show that D(t) ~ t^1/2, while the square of the bracketed term (16pt)^–1. Thus we obtain that Q(t) ~ t^–1/2.

An asymptotic behavior for F(t) follows if we approximate the sum in (4.14) by an integral. Assuming that F(t) ~ t–, we must impose that the integral

does not depend on t. As this is verified when = 1/2, inserting this value into (4.11) leads finally to a = 1/2, the asymptotic value obtained from the numerical evaluation and also found for the DR model [4].

VI Conclusions

In this work we presented an analytical investigation of an abelian SOC model that describes some aspects of rainfall. This model shares many properties of the abelian DR model, but as it allows for the presence of holes within an avalanche cluster, it has some properties of its own.

We were able to obtain an analytical solution for the probability P(s₁; s₀) that a site s₁ topples when one grain is added in any other site s₀ of the lattice. Knowing this probability function is essential for the derivation of many other properties of the model. So we could evaluate the exponent that describes the dependence of the statistics of events as function of event magnitude. Although the value of coincides with that of the DR model in the limit of an infinite lattice, we have shown that for finite size lattice this exponent becomes larger than that for the DR model.

The particular solution for the rainfall model described in this work also opens the probability to analyze the behavior of the presence of holes inside an avalanche cluster. This analysis is currently on the way and will be published elsewhere.

Acknowledgements

The author is much indebted to Dr. S.T.R. Pinho and to Prof. D. Dhar for many stimulating discussions and suggestions. This work has been partially supported by CNPq.

[12] In a private communication, D. Dhar mentioned finding a solution for the particular value e = u = w = 1, that does not fall into the ansatz (3.5).

Received on 1 April, 2003

• [1] D. Dhar, Phys. Rev. Lett. 64, 1613 (1990).
• [2] P Bak, C. Tang, K. Wiesenfeld, Phys. Rev. Lett 59, 381 (1987).
• [3] P. Bak, How nature works: The Science of Self-Organized Criticality, (Copernicus, New York, 1996).
• [4] D. Dhar and R. Ramaswamy, Phys. Rev. Lett. 63, 1659 (1989).
• [5] S. T. R. Pinho and R. F. S. Andrade, Physica A 255, 483 (1998).
• [6] A. Vazquez, e-print cond-mat/0003420 (2000).
• [7] M. Kloster, S. Maslov, C. Tang, Phys. Rev. E 63, 26111 (2001).
• [8] R. Pastor-Satorras and A. Vespignani, J. Phys. A 33, L3, (2000).
• [9] M. Paczuski, K. Bassler, Phys. Rev. E 62, 5347 (2000)
• [10] D. Hughes, M. Paczuski, Phys. Rev. Lett. 88, 054302-1 (2002)
• [11] S. T. R. Pinho, R. F. S. Andrade, A. P. M. Tanajura, and S. C. Fraga, Physica A 314, 405 (2002).

Publication Dates

• Publication in this collection
  11 Nov 2003
• Date of issue
  Sept 2003

History

• Accepted
  01 Apr 2003
• Received
  01 Apr 2003