# Variable predictability in deterministic dissipative sandpile

M. G. Shnirman, A. B. Shapoval

<figure><img src="image_1.png" /></figure>

## To cite this version:

M. G. Shnirman, A. B. Shapoval. Variable predictability in deterministic dissipative sandpile. Non-linear Processes in Geophysics, 2010, 17, pp.85-91. 10.5194/npg-17-85-2010. insu-03605345

HAL Id: insu-03605345

https://insu.hal.science/insu-03605345

Submitted on 11 Mar 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Variable predictability in deterministic dissipative sandpile

M. G. Shnirman<fnref n="1" />,<fnref n="3" /> and A. B. Shapoval<fnref n="2" />,<fnref n="3" />

Received: 24 February 2009 – Revised: 19 January 2010 – Accepted: 11 February 2010 – Published: 24 February 2010

**Abstract.** It is known that some quiescence precedes the strong events in the Bak–Tang–Wiesenfeld sand-pile (Pepke and Carlson, 1994). We introduce dissipation depending on the propagation of the events into this model such that in the constructed model the growth of activity occurs before the strong events. This fact allows the prediction of them in advance with a certain efficiency. This efficiency is variable in time. The best predictability is observed during subcritical time ranges, while the efficiency is definitely worse in the supercritical state.

## 1 Introduction

The attitude to earthquake prediction remains controversial (Wyss, 1997). On the one hand, there exist the prediction algorithms which efficiently forecast strong earthquakes in advance (Kossobokov and Shebalin, 2003). The foreshock activity of middle-size earthquakes underlies these algorithms (Keilis-Borok, 2003). On the other hand, some scientists argue that these algorithms hardly reflect the physics of the seismicity and that the efficiency of the current outcome of the prediction will decline later (Geller et al., 1997).

This discussion is based on the comprehension of the seismic process as a movement of a self-organized critical system of blocks. A typical example of a self-organized critical system is the sandpile introduced by Bak, Tang, and Wiesenfeld (BTW) in Bak et al. (1987). Their model determines the evolution of sand grains on a lattice. The grains are slowly accumulated until their number becomes locally too big. Then they are instantaneously redistributed over the lattice. The redistribution mechanism is conservative inside the lattice and dissipative at the boundary. The slow input and

the quick output balance each other and the system comes to its steady state (Dhar, 1999). Model reviews and open problems can be found at Dhar (2006); Dickman et al. (2000).

According to Pepke and Carlson (1994), (a) strong model events have the *anti-activation* scenario; (b) the adapted earthquake precursors predict these events with a very low efficiency. Further investigation (Shapoval and Shnirman, 2004) of the BTW sandpile gives evidence that its biggest events (which rarely happen) are predictable due to precursors that are unobservable in seismicity. The irregularity of the dissipation (examined by De Menech et al., 1998) and the oscillation of the sand-pile height underlie the prediction. In details, when the number of the grains in the lattice is small the system stays in its subcritical state characterized by a weak dissipation and a rare occurrence of the middle-size events. Then the lattice accumulates “extra” grains and the system comes to the supercritical state. It returns to the subcritical state when a characteristic event happens (De Menech et al., 1998). Just the latter event is predictable due to the quiescence preceding it (Shapoval and Shnirman, 2004). Still, the absence of activation prior to strong events contradicts either the critical self-organization of the seismic process or the predictability of earthquakes based on a certain activation. The “wrong” scenario of strong model events probably happens due to the conservative redistribution of the grains inside the lattice, whereas the fault interactions are dissipative.

We introduce an internal (local) dissipation during the redistribution of the grains into the BTW sand-pile with a central seeding (Wiesenfeld et al., 1990). Dissipation “accompanying” earthquakes seems to be proportional to the volume of the source. The variety of the sand-pile models with dissipative redistribution of the grains exhibits dissipation that is close to the area of the source of model earthquakes. A proper definition of *local* dissipation that ensures the required global features of dissipation (i.e. proportionality to the volume of the source) is not evident. We use the simplest local

<figure><img src="image_3.png" /></figure>

Correspondence to: A. B. Shapoval
(shapoval@mccme.ru)

<footnote n="1">¹Institut de Physique du Globe de Paris, UMR 7154, CNRS, France</footnote>

<footnote n="2">²Finance Academy under the Government of the Russian Federation, Russia</footnote>

<footnote n="3">³International Institute of Earthquake Prediction Theory and Mathematical Geophysics, Russia</footnote>
rule such that dissipation of the full-scale events is unproportionally bigger than dissipation of the small events.

This dissipation leads to the activation scenario of strong events. The constructed model remains predictable but the growth of activity underlies the prediction of strong events. The prediction algorithm corresponds to that for real seismicity (Keilis-Borok, 2003). The highest efficiency of prediction is attained when the system comes to the subcritical state.

# 2 Model

## 2.1 Dynamics

Let $\{(i, j)\}_{i,j=1}^L$ be a square lattice. The set of the integers $\Xi = \{h_{ij}\}_{i,j=1}^L$ is called a “configuration”. These integers are interpreted as the heights of the sand grains on the cell $(i, j)$.

The cell $(i, j)$ is stable if $h_{ij} < H$, where $H = 4$ is the cell instability threshold. If $h_{ij} \ge H$, then the cell $(i, j)$ is unstable. Further, the configuration $\Xi$ is called “stable” if all its cells are stable. Otherwise, the configuration is unstable.

Let a random variable take the values 0, 1, 2, 3 with the probability of $1/4$. Then its $L^2$ independent observations determine the initial configuration $\Xi_0$. Evidently, $\Xi_0$ is stable.

Now we define the mechanism transforming the configuration $\Xi(t)$ appearing at the time step $t$ to $\Xi(t+1)$. Initially, the mechanism adds a new grain onto the central cell $(i_0, j_0)$:

$$h_{i_0 j_0} \rightarrow h_{i_0 j_0} + 1.$$

If the cell $(i_0, j_0)$ remains stable nothing more occurs at this time step $t$. Then the configuration $\Xi(t+1)$ is obtained. If $h_{i_0 j_0} \ge 4$ then sand is redistributed.

We start the definition of the redistribution with one act initiated by any unstable cell $(i, j)$. Let $h_{ij} \ge 4$. Then the unstable cell passes four grains one-by-one to its four nearest neighbours. Several grains can dissipate during this pass:

$$h_{ij} \rightarrow h_{ij} - 4 - D_{ij}, \qquad (1)$$

$$h_{c(i,j)} \rightarrow h_{c(i,j)} + 1 \quad \forall c(i,j), \qquad (2)$$

where $c(i, j)$ is any cell such that it has a common side with $(i, j)$ and $D_{ij}$ is some non-negative integer. Then $D_{ij}$ grains dissipate if the cell $(i, j)$ does not belong to the lattice boundary. The redistribution on the boundary results in one (or two for the corners) additional dissipated grain since the boundary cells have less than four neighbours. It worth reminding that $D_{ij} = 0$ in the BTW sand-pile.

Clearly, the act of the redistribution can lead to the formation of new unstable cells. The redistribution starts in the central cell $(i_0, j_0)$ and continues until a stable configuration occurs (during the redistribution one cell can become unstable several times). The final configuration is just $\Xi(t+1)$. Then the next time step begins.

Two time scales are defined in the model. The grain addition occurs during slow time (the word “slow” is usually omitted). The redistribution is associated with quick time.

## 2.2 Dissipation

It remains to define $D_{ij}$. Let $z_{ij}$ be the counters of the cell instability. It is supposed that $z_{ij} = 0 \forall i, j$ at the beginning of any time step. Once the cell $(i, j)$ becomes unstable the counter $z_{ij}$ increases by one. The threshold value $z^*$ is fixed for all the counters $z_{ij}$. By definition, put

$$D_{ij} = \begin{cases} 0, & \text{if } z_{ij} < z^*; \\ d^*, & \text{if } z_{ij} \ge z^*, \end{cases} \qquad (3)$$

where $d^*$ is some natural number. Thus in addition to a boundary dissipation we define an internal dissipation depending on the redistribution of the grains.

Our model depends on two parameters: $z^*$ and $d^*$. The value $z^*$ determines when the dissipation is switched on. The value $d^*$ stays for the number of the dissipating grains during one act of the redistribution.

The value of the instability threshold $H$ does not influence the model dynamics. We fix $H = 4$ to make easier the comparison with the BTW sand-pile. In this case the heights can become negative due to dissipation. However they are still bounded from below. Whence, moving $H$ higher makes all the heights positive. Interpreting $h_{ij}$ as the local stress we keep in mind the model with sufficiently big $H$.

If the redistribution occurs at some time step then this process is called an event. Its size is the number of the unstable cells appeared during the redistribution and counted with regard for multiplicity.

The model differs from the BTW sand-pile in two ways. First, new grains are added not in cells chosen at random but in the lattice center only (as discussed in Wiesenfeld et al., 1990). This makes the dynamics deterministic. Secondly, the dissipation $D_{ij}$ differs from zero.

# 3 Power recurrence law

The BTW sand-pile has gained popularity due to its power recurrence law. We draw the corresponding plot for the constructed model. Let a catalogue be the set of the consequent events. The catalogue has to be sufficiently extent and remote from the initial configuration such that the system is able to attain the steady state. Suppose $\Delta s$ is some fixed number being insignificantly bigger than one. Let $F(s)$ be the number of the catalogued events, whose size lies in $\sigma \in [s/\Delta s, s\Delta s)$, divided by the extent of the catalogue (in other words, by the number of the grains added). If $F(s)$ is a power function then its graph is linear in the log-log scale.

Figure 1 introduces the function $F(s)$ for the different values of the parameters $z^*$ and $d^*$. The typical pattern for the function $F(s)$ in the log-log scale consists of the almost linear part followed by the bump and the abrupt break down. To check whether this pattern or at least its linear part is conserved as $L \to \infty$ one has to simulate the models for several bigger $L$s. In the paper we leave this problem aside. The big
<figure><img src="image_4.png" /><figcaption>**Fig. 1.** The frequency $F(s)$ vs. the size $s$ for the different models; "$z^* = infinity$" stays for the BTW model with a central seeding; the slopes of the dashed lines are -0.12 and -0.15; the catalogues' extent is $5 \times 10^5$; $L = 128$; $\Delta s = 10^{0.1}$.</figcaption></figure>

values of the parameter $z^*$ mean that the dissipation $D_{ij}$ is (almost always) equal to zero (according to Eq. 3). Then the model dynamics agrees with the BTW sand-pile with grain addition in the central cell. If $z^*$ is small the dissipation is switched on so early that the strong events are hardly realized. Therefore the graph for $z^* = 4$, $d^* = 1$ lies significantly lower than the others.

Note that the function $F(s)$ exhibits locally cumulative size distribution. This increases the slope of the power part of the graphs by 1 (in the log-log scale) comparatively with the direct size-frequency plot given, for example in Wiesenfeld et al. (1990) for the BTW model with the central seeding.

# 4 Prediction

## 4.1 Precursor “Power Size” (PSi)

There are prediction algorithms forecasting strong earthquakes in advance (Keilis-Borok and Kossobokov, 1990; Keilis-Borok and Rotwain, 1990; Kossobokov and Shebalin, 2003; Shebalin, 2006). The precursors of the strong earthquakes underlying these algorithms quantitatively describe the increase of the middle-size earthquakes prior to the strong earthquakes. We are going to adapt these precursors to the big model events. The event is called “target” if its size is bigger than some $s_0$. The target events correspond to the strong earthquakes. Suppose [$s_-, s_+$], $s_- < s_+ < s_0$, is the size interval of the middle-scale events, $w$ is the length of the

<table><caption>**Table 1.** Optimal values of the parameters fixed for different pairs ($d^*, z^*$).</caption><thead><tr><th>d*</th><th>z*</th><th>lg(s<sub>-</sub>)</th><th>lg(s<sub>+</sub>)</th><th>α</th><th>w</th><th>T</th></tr></thead><tbody><tr><td>1</td><td>4</td><td>2.80</td><td>3.80</td><td>0.25</td><td>300</td><td>100</td></tr><tr><td>1</td><td>5</td><td>3.00</td><td>4.00</td><td>0.25</td><td>200</td><td>100</td></tr><tr><td>1</td><td>9</td><td>3.45</td><td>4.45</td><td>0.25</td><td>70</td><td>300</td></tr><tr><td>1</td><td>10</td><td>3.50</td><td>4.50</td><td>0.25</td><td>100</td><td>100</td></tr><tr><td>1</td><td>20</td><td>3.80</td><td>4.80</td><td>0.25</td><td>100</td><td>300</td></tr><tr><td>10</td><td>10</td><td>3.25</td><td>4.25</td><td>0.75</td><td>350</td><td>100</td></tr></tbody></table>

sliding window, and $s(t) \in [s_-, s_+]$ is the size of the middle-scale event occurred at the time step $t$. Then the functional

$$ \Psi_{\alpha}(t) = \sum_{k=t-w}^{t-1} (s(k))^{\alpha}, \qquad (4) $$

where $\alpha$ is an appropriate power, is a precursor (of the target events) measuring the occurrence of the middle-scale events.

We name this precursor “Power Size” (PSi) as well as the prediction algorithm based on this precursor. Prediction algorithm PSi includes the calculation of the precursor and a rule switching *alarms* on and off. As soon as $\Psi_{\alpha}(t) > \Psi^*$ (for appropriate fixed $\Psi^*$) the algorithm expects a target event to occur during the next $T$ time steps. More precisely, let $t_{on}$ be any time step specified by $\Psi_{\alpha}(t_{on}) > \Psi^*$. Then $t_{off} = t_{on} + T$ is assigned to $t_{on}$ if the target events are absent on $[t_{on}, t_{on} + T]$. Otherwise by $t_{off}$ denote the step of the first target events. Then the union of all $[t_{on}, t_{off}]$ (which possibly intersect one another) forms a *time of increased probability* (TIP) of the target events. The target events are said to be *predicted* if they occur during TIP.

## 4.2 Prediction efficiency

The ratio of the unpredicted events and the TIP extent naturally describe the prediction efficiency (Keilis-Borok, 2003; Molchan, 2003). Suppose $n$ is the ratio of the unpredicted events, $\tau$ is the ratio of TIP (in other words, $\tau$ is the total TIP extent divided by the catalogue length), and $\varepsilon$ is $n + \tau$. Then $\varepsilon$ is called the *loss* of the algorithm. According to Molchan (2003), the loss $\varepsilon$ is close to 1 for a random prediction. The *efficiency* of the prediction is $1 - \varepsilon$.

## 4.3 Choice of parameters

The parameters $s_-$, $s_+$, $\alpha$, $w$, and $T$ and $\Psi^*$ are adjusted on some “learning catalogue” to minimize $\varepsilon$. In fact, $\varepsilon$ depends on the parameters in a complex way. Therefore a parallelepiped lying in the parameter space is detected “by hands”. This parallelepiped is gridded and each node (which is a point in the parameter space) is examined. The node generating the least $\varepsilon$ gives the values of the parameters that are called optimal and fixed (Table 1). Then the algorithm is applied to another catalogue with all these values fixed.
<figure><img src="image_5.png" /><figcaption>Fig. 2. The loss $\epsilon$ with error bars vs. the threshold $z^*$ switching the dissipation on for $d^* = 1$. The numbers are the values of $\log(s_0)$ (the lowest size of the target events).</figcaption></figure>

The low boundary $s_0$ of the target events is not adjusted. Aiming at the best efficiency we choose $s_0$ as big as possible (this idea works for the BTW sand-pile; Shapoval and Shnirman, 2004). The fixed values of $s_0$ (being written in Fig. 2 for $L=128$) ensure approximately the same frequency of the target events in the catalogues for different $z^*$. Despite these values of $s_0$ are bigger than the sizes shown in Fig. 1 the number of the target events in the catalogue is big enough to get statistically significant results.

## 4.4 Prediction outcome

Figure 2 introduces the loss $\epsilon$ for the prediction algorithm corresponding to the dissipation $d^* = 1$ and several values of $z^*$. If $z^*$ is too big algorithm PSi does not predict the target events. This agrees with the BTW sand-pile. Small values of $z^*$ ($z^* = 4$ in Fig. 2) lead to a weak predictability. The strong events are extremely rare in such models.

The best predictability is observed for intermediate $z^*$ (Fig. 2). These values of $\epsilon$ are far from 1. This gives evidence that the prediction results are not random. Moreover, the efficiency increases whenever the dissipation goes up. We do not support this statement for all $z^*$ but give the detailed analysis of the predictability in the model determined by $z^* = 10$, $d^* = 10$.

So, fix $z^* = 10$, $d^* = 10$. Algorithm PSi is applied to the catalogue sampled during $[2.5 \times 10^6, 7.5 \times 10^6]$ time steps ($L=128$; $5 \times 10^6$ grains are added on the lattice). It keeps 2601 target events. The algorithm predicts 2185 events while the alarm continues about a third of the catalogue's extent. In other words, the outcome of the prediction is $n = 0.16$, $\tau = 0.35$, $\epsilon = 0.51$. To check this result other catalogues of the same length are generated for different initial configurations. The values of $\epsilon = n + \tau \approx 0.5$ are conserved.

The parameters of the algorithm PSi influence the efficiency in a different way. The most principal parameters are $w$ and $T$. On the contrary, the parameter $\alpha$ determining the power in the functional $\Psi_\alpha$ defined in Eq. (4) weakly influences the efficiency. We claim that the loss of the algorithm PSi is less than 0.54 as $\alpha \in [0.1, 1]$.

## 4.5 Role of dissipation

We have determined the nonlinear dissipation depending on the propagation of the model events. This dissipation weakens the middle-scale events preventing their propagation. Therefore it takes the series of the middle-scale events to transport the grains to the boundary. Only then the full-scale event happens. Hence it can be predicted. This scheme gives the explanation of the reported prediction efficiency.

## 4.6 Efficiency variability

By definition, put

$$h(t) = L^{-2} \sum_{i,j=1}^{L} h_{ij}(t),$$

where the values of $h_{ij}(t)$ are taken at the end of the time step $t$. Further, suppose $\langle h \rangle(t)$ is the mean of $h(t)$ over previous $N_s=50000$ time steps and $\langle \epsilon \rangle(t)$ is the loss $\epsilon$ calculated on $[t - N_s, t]$. Then the oscillations of $\langle \epsilon \rangle(t)$ are rather big (Fig. 3). The loss can be doubled (from 0.35 to 0.70) due to the choice of the prediction interval. Is this variability connected with the sandpile height, which changes significantly too (Fig. 3)?

By $\rho(t)$ we denote the correlation function of $\langle \epsilon \rangle(t)$ and $\langle h \rangle(t)$. It is calculated on the intervals $[t - N_b, t]$, $N_b=250000$. Numbers $N_s$ and $N_b$ balance two opposite requirements. On the one hand, the used intervals have to be small such that the model properties described in terms of $h(t)$ do not change a lot. On the other hand, if the intervals are too small then their number of the target event remains insufficient for meaningful conclusions about the prediction efficiency.

According to Fig. 3, lower panel, there exist extremely long time ranges (hundreds of thousand time steps) exhibiting the correlation $\rho$ of the constant sign. Still the mean of the correlation function is close to zero.

## 4.7 Sub- and supercriticality

Our prediction method is expected to be appropriate near the steady state (Sornette, 2002). The steady state consists of a large number of configurations with approximately constant height of the sand pile. This height being close to 1.88 in the model investigated ($d^* = 10$, $z^* = 10$) can be interpreted as a critical height. Nevertheless now and then $\langle h \rangle(t)$ becomes too small (Fig. 3) since the biggest events are strongly dissipative. Hence the height goes far away from its critical level.
<figure><img src="image_6.png" /><figcaption>**Fig. 3.** A typical evolution of the smoothed height's mean, prediction efficiency (upper panel), and their correlation function $\rho(t)$ (lower panel) for $d^* = 10$, $z^* = 10$; $\langle h \rangle(t)$ and $\langle \varepsilon \rangle(t)$ are the average of $h(t)$ and, respectively, the loss of the prediction, over $[t - 5 \times 10^4, t)$. The dotted curve stays for $\rho(t)$ on the cutting off intervals with too low sand level. The solid curves stay for the other intervals. The mean of the solid curves (of $\langle h \rangle(t)$) is approximately 0.27 (1.88, dashed lines). Zero on the time axes corresponds to the beginning of the catalogue.</figcaption></figure>

We want to analyze only the parts of the correlation function in which the average height is greater than or equal to an arbitrary threshold, $\langle h \rangle(t) < h^*$, because a part of our hypothesis is that better prediction is not possible when the average height is low (corresponding to few events). Consequently, for each time step $t'$ where $\langle h \rangle(t') < h^*$, we identify the corresponding time interval (in the past) over which $\langle h \rangle(t')$ was previously defined, namely $[t' - N_s, t']$. We want to exclude all the points that lie in this interval from the calculation of any future values of the correlation function $\rho(t)$. First, notice that time steps within this “forbidden” interval will be used in the calculation of $\langle h \rangle(t)$ for all values of $t$ in $[t', t' + N_s]$. Second, notice that time steps within the interval $[t', t' + N_s + N_b]$ will determine $\rho(t)$ using values of $\langle h \rangle(t)$ which themselves use data from the “forbidden” interval. Therefore, for any value of $t'$ for which $\langle h \rangle(t') < h^*$, we mark the interval $[t', t' + N_s + N_b]$ with a dotted line in the correlation plot of Fig. 3 to signify regions with poor predictability due to low system average height (or low system mass). The remaining part of the correlation plot is marked with thick solid line to identify regions with better predictability.

According to Fig. 3, the remaining part of $\rho(t)$ (plotted by the solid line) is more positive than negative. The mean $\bar{\rho}$ of this $\rho(t)$-part shown in the figure is 0.27 (for reliable conclusion the mean is calculated for two other catalogues; the values are 0.25 and 0.30). Hence the height and the efficiency fluctuate co-directionally if the height is close to its critical level. Then the growth of the height usually increases the loss $\varepsilon$ of the prediction algorithm. If the height is near its critical level then the height's growth pushes the system to the supercritical state. Thus the co-directional oscillations of $\langle \varepsilon \rangle(t)$ and $\langle h \rangle(t)$ give an implicit evidence that the predictability is worse in the supercritical state.

The value of $h^*$ has to be sufficiently big to study the system near its critical level of height. However excessively big $h^*$s eliminate the graph of $\rho(t)$ completely. It is fixed $h^* = 1.7$ whenever the mean of $h(t)$ is approximately 1.88. (We check the values of $h^*$ for three considered catalogues; $h^* = 1.65$ implies $\bar{\rho} = 0.14; 0.25; 0.30$ and $h^* = 1.75$ does $\bar{\rho} = 0.26; 0.27; 0.21$; this calculation gives evidence of a certain stability of the results with respect to $h^*$).

The direct calculation of $\varepsilon$ in the subcritical and supercritical state leads to the values 0.44 and 0.57 respectively. The
calculation corresponds to Fig. 3’s all the intervals such that $h(t)$ belongs to [1.7, 1.85] and [1.92, 2.0] during at least 500 consequent time steps.

The following idea can explain the observed variability. If the number of the grains on the lattice is too big then even minor changes in the configurations are able to generate a strong event. A long formation process is not necessary. A strong event can occur “at any moment”. Therefore the predictability is weaker in the supercritical state. Furthermore, extremely low pile height assures that the growth of activity in a lattice part does not lead to a strong event because the other part lacks the sand grains. A strong event would occur later (and remain unpredictable) when the lattice accumulates some additional grains and a new wave of the redistribution attains the overloaded lattice part. This leads to a low prediction efficiency. Finally, whenever the height comes to its critical level from below there exists the formation process of the strong events involving the increase of activity.

## 4.8 Accuracy of results

The properties of the prediction developed in the paper can be sensitive to the choice of the parameters. The applied scheme searches the minimum of $\varepsilon$ through a limited number of points lying in the many-dimensional parameter space. On the one hand, the loss $\varepsilon$ of the prediction as the function in any particular parameter has a V-shape on the grid if the values of the other parameters are fixed as reported in Table 1 (we do not accompany this fact by figures). If $\varepsilon$ depends on the parameters sufficiently smoothly then the obtained node of the parameter grid is close to the global minimum of $\varepsilon$. On the other hand, $\varepsilon$ as a complex non-linear multivariable function can have irregularities that are not described by the values calculated on the grid's nodes. Whence the global minimum of $\varepsilon$ can be passed through. If the latter is true (which is unlikely to happen from our point of view) then the $(n, \tau)$-outcome of the prediction can be only more efficient than that found in this research. However the other prognostic properties of the model (in particular, the variability of the prediction) has to be verified for the global minimum of $\varepsilon$. Anyway, these properties are valid given a natural construction involving the node-by-node examination of a reasonable parameter grid.

## 5 Conclusions

We introduce a non-linear dissipation of the sand grains during their redistribution into the BTW sand-pile with a central seeding. Several features of seismicity are built in the model. They are two time scales and the locality of the mechanism running the propagation of the events. The model dynamics follows other seismic feature:

- the power recurrence law,

- the predictability of the strong events based on the activation,

- non-stationarity of the prediction.

Two last properties (assured by the introduced dissipation) separate the constructed model from the BTW sandpile and its simple modifications, where a certain *quiescence* precedes the strong events (Pepke and Carlson, 1994; Shapoval and Shnirman, 2009). Thus the constructed model realizes the predictability of the seismic process based on the *activation* inside the class of the self-organized critical systems. Earthquakes can be preceded by a certain combination of activation and quiescence (Huang et al., 1997) but the model construction of this phenomenon is the subject of a separate study.

In the developed model the dissipation plays a central role. The seismic process is characterized by the number of failures being inversely proportional to the failure area, while the dissipation is proportional to the failure volume. That is why the strongest earthquakes are accompanied by the most noticeable dissipation of the seismic process. Hence a model analogue of the real dissipation has to be a function on the event's size, which grows more quickly than linearly. The best choice of this function has not been discovered yet. Therefore we focus on the simplest nonlinear function emphasizing the dominating dissipation of the strongest events.

The model system oscillates in its steady state such that the typical time range of the oscillations is extremely big. The predictability of the strong events is definitely better in the subcritical state than in the supercritical state. In the model terms, a local non-stationarity of the seismic process can result in a temporal break down of the real prediction.

*Acknowledgement.* We are thankful to the anonymous reviewers for their valuable suggestions and comments. We acknowledge a partial support from Russian Foundation for Basic Research (Grants No. 08-05-00215-a, No. 10-06-00282-a, and No. 08-01-00784-a). This work was supported by the Institut de Physique du Globe de Paris (IPGP contribution 2479).

Edited by: U. Feudel

Reviewed by: two anonymous referees

<figure><img src="image_7.png" /></figure>

The publication of this article is financed by CNRS-INSU.
References

Bak, P., Tang, C., and Wiesenfeld, K.: Self-Organized Criticality: An Explanation of 1/f Noise, Phys. Rev. Lett., 59, 381–384, 1987.

De Menech, M., Stella, A. L., and Tebaldi, C.: Rare Events and Breakdown of Simple Scaling in the Abelian Sandpile, Phys. Rev. E, 58, R2677–R2680, 1998.

Dhar, D.: The Abelian Sandpile and Related Models, Physica A, 263, 4–25, 1999.

Dhar, D.: Theoretical Studies of Self-Organized Criticality, Physica A, 369, 29–70, 2006.

Dickman, R., Munoz, M., Vespignani, A., and Zapperi, S.: Paths to Self-Organized Criticality, cond-mat/9910454v2”, 2000.

Geller, R. J., Jackson, D. D., Kagan, Y. Y., and Mulargia, F.: Earth-
quakes cannot be predicted, Science, 275, 1616–1617, 1997.

Huang, Q., Sobolev, G. A., and Nagao, T.: Characteristics of the seismic quiescence and activation patterns before the $M = 7.2$ Kobe earthquake, January 17, 1995, Tectonophysics, 337, 99–116, 2001.

Keilis-Borok, V. I.: Fundamentals of Earthquake Prediction: Four Paradigms, in: Nonlinear Dynamics of the Lithosphere and Earthquake Prediction, edited by: Keilis-Borok, V. I. and Soloviev, A. A., Springer-Verlag, Heidelberg, pp. 1–36, 2003.

Keilis-Borok, V. I. and Kossobokov, V. G.: Preliminary activation of seismic flow: Algorithm M8, Phys. Earth Planet. Int., 61, 73–83, 1990.

Keilis-Borok, V. I. and Rotwain, I. M.: Diagnosis of time of increased probability of strong earthquakes in different regions of the world: algorithm CN, Earth Plannet. Inter., 61, 57–72, 1990.

Kossobokov, V. and Shebalin, P.: Earthquake Prediction, in: Non-linear Dynamics of the Lithosphere and Earthquake Prediction, edited by: Keilis-Borok, V. I. and Soloviev, A. A., Springer-Verlag, pp. 141–208, 2003.

Lübeck, S., Rajewsky, N., and Wolf, D.E.: A deterministic sandpile automaton revisited, Eur. Phys. J. B., 13, 715–721, 2000.

Molchan, G. M.: Earthquake Prediction Strategies: A Theoretical Analysis, in: Nonlinear Dynamics of the Lithosphere and Earthquake Prediction, edited by: Keilis-Borok, V. I. and Soloviev, A. A., Springer-Verlag, pp. 209–238, 2003.

Pepke, S. L. and Carlson, J. M.: Predictability of Self-Organizing Systems, Phys. Rev. E, 50, 236–242, 1994.

Shapoval, A. B. and Shnirman, M. G.: Strong events in the sand-pile model, Int. J. Mod. Phys. C, 15, 279–288, 2004.

Shapoval, A. B. and Shnirman, M. G.: Earthquake Precursors used for predicting the largest events in the avalanche formation model, Izvestiya, Phys. Solid Earth, 45, 406–413, 2009.

Shebalin, P.: Increased correlation range of seismicity before large events manifested by earthquake chains, Tectonophysics, 424, 335–349, 2006.

Sornette, D.: Predictability of catastrophic events: material rupture, earthquakes, turbulence, financial crashes and human birth, Proc. Nat. Acad. Sci. USA, 99, 2522–2529, 2002.

Wiesenfeld, K., Theiler, J., and McNamara, B.: Self-organized criticality in a deterministic automaton, Phys. Rev. Lett., 65, 949–952, 1990.

Wyss, M.: Cannot Earthquakes Be Predicted? Science, 278, 487–490, 1997.