1 Introduction

Self organized critical (SOC) systems [1, 2] have been extensively studied over the past forty years. The criticality of such systems is not based on fine tuning of the system parameters. These systems are slowly driven, and in the critical state result in dissipation, the magnitude of which is power law distributed and scale free. SOC models, therefore, have been appropriate systems for the study of extreme events. An important question which arises is with regard to what degree these extreme events in SOC models can be predicted using past data.

Natural systems such as earthquakes [1, 3, 4], wildfires [5], solar flares [1] show long-range spatial correlations and power-law behavior. These characteristic properties are well reproduced in simple SOC lattice models; this has led to a new paradigm in the prediction of extreme events in natural systems [6]. Of the SOC lattice models, sandpile models are the most convenient to deal with [7]. Multiple studies have used sandpiles and other SOC models to discuss the issue of predictability with regard to natural systems [8,9,10]. In [8], a nonconservative SOC lattice model is shown to be equivalent to the well established Burridge Knockoff model for faults in earthquakes [11, 12]. In [9], the authors identify precursors to large events in a variety of SOC models, analogous to seismic signatures before large earthquakes; and use these to plot prediction success rate curves. In [10], a new sandpile-like automaton model is introduced which successfully reproduces the Gutenberg Richter law for energy distribution in earthquakes, as well as other empirical data of earthquakes from seismogenic regions.

The notion of predictability in sandpiles varies amongst different works. In [13], the sandpile is said to be predictable if, given any initial state of the sandpile, there exists an algorithm to determine the final recurrent state, such that the algorithm is faster than direct time evolution of the system. In this sense, they conjecture that the sandpile model is not predictable. A different kind of problem is short term prediction. Given a past time series of occurence of extreme events, one tries to predict whether an extreme event is likely in the next few timesteps [14, 15]. We will use the term ‘predictability’ in this second sense. Note that a time-evolving process can often be decomposed into a deterministic part, which is perfectly predictable, and a random part wherein possible outcomes are probabilistically weighted. A ‘blind’ prediction algorithm which randomly guesses outcomes, predicts correctly a fraction of the events. We state that an event is predictable if there exists a strategy which can can correctly predict a larger fraction of outcomes, than such a blind algorithm.

Next timestep prediction is based on a decision variable constructed out of present or past values of observables. In [14], the decision variable is the fraction of almost-unstable sites on the sandpile. In [16], the decision variable is constructed with a weighted sum of past avalanche sizes, where the kernel is optimized for prediction efficiency. In both studies, they find that in the thermodynamic limit, sufficiently ‘large’ events (to be defined later) in Manna sandpiles [17] can be predicted with a good success rate. However, BTW sandpiles do not seem to be predictable. A similar strategy is used in [15], where subsequent predictions for extreme events are decreased if they occured recently in the past. This predictability, based on anticorrelation of extreme events, is based on a finite size effect, and vanishes in the thermodynamic limit. Both these strategies depend on continuous internal parameters. As a ‘measure’ of the quality of prediction, predictability curves are plotted for a wide range of these parameters.

We will discuss here the prototypical Directed sandpile model [7, 18], because it is more tractable than the BTW or Manna sandpiles. This allows us to analytically calculate predictability curves (defined in detail later) and scaling forms in many cases. One can make many different kinds of ‘predictions’ [19, 20], including predictions about the strength of the event, or given past data, provide the full a posteriori probability distribution of event sizes. We consider the simplest setting, where the prediction is a binary variable, in the form of an alarm, which is turned on if the extreme event is likely, and off otherwise. Of course, the strategy to turn the alarm on or off at a time (t+1), takes the past data for times (1, .., t) as input. We compare the predictability of two different strategies. The first is a waiting time strategy, which exploits the temporal anti-correlation of large avalanches. For an (L \times M) directed sandpile, with (M \sim L^{1/2}), we find that large avalanches are predictable to some extent. More importantly, contrary to [15, 16], this predictability is not a finite size effect and persists in the thermodynamic limit. Our second strategy makes use of measuring the excess density in the model, similar to [14]. We find that with such a strategy, avalanches are somewhat predictable for finite L, but predictability decreases with L and vanishes in the thermodynamic limit.

This paper is structured as follows. In Sec. 2, we introduce the directed sandpile model, and general terminology related to the prediction problem. In Sec. 3 we first discuss the efficiency of the waiting time strategy to predict large events. In Sec. 3.1 we show, using MC simulations that predicting rare events everywhere globally in the lattice is not possible under this strategy. In Sec. 3.2 we obtain the waiting time distributions for events occurring in a local detection window. We present analytic expressions for the average waiting times. We plot predictability curves (defined later) using Monte Carlo data, and theoretically find expressions for the prediction efficiency as a function of event size and system size, in the large L limit. We show that there is a finite non zero gain in prediction efficiency in the thermodynamic limit. Next, we consider a different kind of prediction strategy based on the excess density in the sandpile. In Sec. 4.1 we use the known scaling properties of large avalanche events to exactly determine the predictability curves for the prediction strategy. We explicitly show that predictability diminishes with lattice size, and vanishes in the thermodynamic limit, and verify this using MC simulations. In Sec. 4.2, we first derive the conditional probabilities for size events, given an excess density, and use this to get exact expressions for the prediction curves. Again, we show that some non zero prediction advantage exists only for finite lattice size. This is also verified by MC simulations. Finally, in Sec. 4.3, we compare these results with previously established numerical studies of other sandpile models.

2 Model

We consider the directed sandpile model [18] on a diagonal square lattice. The vertical height is L, and there are periodic boundary conditions in the horizontal direction. Each horizontal cross section has M sites (Fig. 1). At each lattice site (\vec {r}) there is an associated number of sand grains (z_{\vec {r}} \ge 0). If (z_{\vec {r}} \ge 2), the site (\vec {r}) is said to be unstable.

Fig. 1

An (L\times M) directed sandpile. Grains are added in the top layer. Unstable sites topple and grains move to the downward diagonal neighbours. Grains flow out from the bottom layer

At each time step, one sand grain is added at random to one of the sites (\vec {r}_0) in the top layer, so that (z_{\vec {r}_0} \rightarrow z_{\vec {r}_0} +1). If this makes the site unstable, the site is ‘toppled’, and 2 sand grains are removed from (\vec {r}_0), and 1 grain is transferred to each of the two immediate downward diagonal neighbours. In this manner, all unstable sites are successively toppled, until the sandpile reaches a stable configuration. Note that the toppling of a site in the bottom row results in 2 grains flowing out of the sandpile. The toppling moves are Abelian. A series of topplings of sites one below another, in the same timestep, is called an avalanche. The tractability of this model arises from the fact that each site topples at most once in an avalanche, and the avalanche progresses in a directed manner to the bottom of the sandpile.

Our analysis involves next timestep prediction of ‘large’ avalanches in directed sandpiles. There are different ways to quantify the ‘size’ of an avalanche. We denote the total number of toppled sites in an avalanche as (S_t) and the vertical height of the avalanche as (h_t). The number of topplings in the bottom row is called as the drop number, (f_t). The net outflow of particles in a timestep t is (2f_t). We are given a time series of past data about sandpile events, upto time T. This could include observables (such as excess density (\epsilon _t)) or sizes of events. Given this data we want to predict an extreme event at time (T+1). We make binary predictions, wherein we turn on an ‘alarm’ if we expect an extreme event to happen with appreciable probability, else we turn it off.

Notation : (\left( E_1, E_2, ..,E_T \right) ) is a binary time series where (E_i = 1) if there is a large event at time i, and (E_i = 0) otherwise. Similarly, (\left( X_1, X_2, ..,X_T \right) ) is the binary series of past predictions. At some time t, we set (X_{t+1} = 0) (alarm off) or (X_{t+1} = 1) (alarm on). Four possible outcomes arise with regard to our prediction : it is a true positive (tp), a true negative (tn), a false positive (fp) or a false negative (fn). The fractions for each of these are given by :

$$\begin{aligned} {\text {tp}}&= T^{-1}\sum _{t=1}^{{T}\delta _{E_t, 1}\delta _{X_t, 1} , , , , {\text {tn}} = T}{-1}\sum _{t=1}^{{T}\delta _{E_t, 0}\delta _{X_t, 0} , , , , {\text {fp}} = T}{-1}\sum _{t=1}^{{T}\delta _{E_t, 0}\delta _{X_t, 1} , , \nonumber \ {\text {fn}}&= T}{-1}\sum _{t=1}^{T}\delta _{E_t, 1}\delta _{X_t, 0} , . \end{aligned}$$

(1)

The prediction variables ({ X_t}) are determined by the prediction strategy. Many prediction strategies are dependent on an intrinsic parameter (\sigma ) (a threshold) which can be chosen to optimize the quality of prediction. A popular class of curves that are plotted to judge the quality of predictions about a system are called ROC curves [21]. In such curves, the false positive rate (y = \frac{{\text {fp}}}{{\text {tn +fp}}}) is plotted against the false negative rate (x = \frac{{\text {fn}}}{{\text {tp + fn}}}), for a wide range of the parameter (\sigma ). The resulting curve is expected to be concave and strictly decreasing.

We prefer to use an equivalent, but more direct characterization, where on the vertical axis we plot the fraction of positive events that are correctly predicted (c(\sigma ) := \frac{{\text {tp}}}{{\text {tp+fn}}}). This is plotted against the average fraction of time for which the alarm is on (\nu (\sigma ) := {\text {tp + fp}}). We will refer to this as an SROT (Success Rate versus On-Time) curve. The goal is to get the best predictability for the least alarm on time. For an entirely unpredictable (uncorrelated) system, the SROT curve is a straight line (c = \nu ). For a partially predictable system, the SROT curve is expected to be convex and strictly increasing for (\nu \in [0,1]).

3 Strategy using waiting times

We assume that the signal (time series) containing information about extreme events is a stationary stochastic process for adiabatic driving rates (i.e. we wait for the configuration to stabilize before adding the next sand grain at the top). In the directed sandpile, there is an anti-correlation between extreme events. For instance, a large avalanche at (t = 0) drains out a lot of particles from the sandpile. So, a future avalanche propagating through the same region is unlikely to be as big. Our strategy simply makes use of this low probability window to switch off the alarm. The strategy makes no use of any specific observables, but only the general statistical behaviour of the time series. We assume this statistics, in particular the distribution of waiting time intervals, is known. We define the following kinds of events as large: for toppling numbers, (S \ge S^* = AL^{3/2}), for heights (h \ge h^* = zL ), and for drop numbers (f \ge f^* = AL^{1/2}). Here, A, z are (\mathcal {O}(1)) parameters. We call these events as ‘large’ events because such avalanche sizes are much larger than the average values, which scale as (\overline{S} \sim L, \overline{h} \sim L^{1/2}) and (\overline{f} \sim L^{0}) for the directed sandpile.

Our strategy is as follows. Fix a time (\tau =\tau ’ L^{1/2}) which is called the off time interval. Suppose there is a large event at time (t_0 = 0). We switch off the alarm for the subsequent (\tau ) timesteps, and switch it back on at time (\tau + 1). The alarm now remains on until another large event occurs.

We will consider the directed sandpile on an (L\times mL^{1/2}) lattice. For large enough m, the local properties of the model are expected to be independent of m. This is because laterally distant sites in the model do not interact. The only effect of changing m is that the effective timescale for local events is proportional to m.

We vary the off time parameter (\tau ), and use simulation data to obtain (c(\tau )) and (\nu (\tau )). We then plot SROT curves.

3.1 Monitoring large events globally

Consider the (m \rightarrow \infty ) limit for our (L \times mL^{1/2}) lattice. The time series available includes sizes of all avalanche occurring throughout the lattice. Henceforth, we call this the ‘global waiting time strategy’, and call the monitored events as ‘global events’. Numerically, it is adequate to take (M = L), because avalanches have widths that scale as (L^{1/2}). We perform Monte Carlo simulations, and generate a time series (length (10^8)) of events for various L.

3.1.1 Waiting time distributions for global events

We plot the normalised waiting time distributions for ‘global events’ of the form (S \ge AL^{3/2}) and (h \ge zL). Denote these distributions as (w_g^{S}(t)) and (w_g^{h}(t)). For a (160 \times 160 ) lattice and the above-mentioned time-series, the plots are in 2a and 2b respectively.

For large A values the distribution is almost flat for a large range of t, because localized events far from the area of interest act as an effective noise for the waiting time distribution. The typical width of a big avalanche is (\mathcal O\left( L^{1/2}\right) ), so a ‘local’ region on the lattice around a point of interest can be defined as an (L\times L^{1/2}) window. On average, the frequency of events outside such a local window is (ML^{-1/2}) times that of events in the window. So, almost all the events in the ‘global’ time-series are uncorrelated in the large M limit, and the global waiting time distribution is simply an exponential [15, 22, 23]. This property is also seen in Fig. 2b. This global strategy cannot be used for prediction, because the flat nature of the global waiting time distribution does not give us a good threshold to switch off the alarm.

We also plot the average waiting time for toppling number events, (\overline{t_g^S}) (as a function of A), and height events, (\overline{t_g^h}) (as a function of z). The plots are in Fig. 3a, 3b, with (L = 90,120, 160). As one would expect, we get a scaling collapse when the average waiting times are scaled by (L^{-1/2}).

Fig. 2

Global event waiting time distributions for different A for (a) Toppling number variables. Note that the x axis is (\ln t), this has been used because the distributions are extremely heavy tailed. (b) Height variables plotted for different thresholds z

Fig. 3

Scaled average waiting times: (a) For toppling number events (b) For height events. Note that (\ln \left( L^{-1/2} ,\overline{t}_g^{h}\right) ) is plotted on the y axis

3.1.2 SROT plots

It is clear that the predictability for global events is not going to be very good. So, instead of making the usual c vs (\nu ) plot, we plot (c-\nu ) vs (\nu ) using the Monte Carlo data. We sample a wide range of waiting times (\tau ), and lattice sizes (L = 90,120,160). The curves are shown in Fig. 4a, 4b. As expected, the prediction is not very effective; the gain in predictability is of (\mathcal O\left( 10^{-2}\right) ) for the considered lattice sizes.

Fig. 4

Global waiting time strategy SROT plots for (M = L)

3.2 Local waiting time strategies

Consider the case where we have a local ‘detector’, and only take into account events that occur within a finite distance of it. Given the noise in data with global events, it is clear that this will lead to better predictions. We use a lattice of size (L \times mL^{1/2}) with large enough m, and periodic boundary conditions. Within this, we consider a smaller vertical window of height L and width (L^{1/2}). We only take into account the outflow from the bottom row in this window. For simplicity, in this section we discuss the case when the prediction is based on drop number events (f \ge AL^{1/2}). Note that only the outflow within the window is detected.

3.2.1 Waiting time distributions for drop number variables

We first analyse the waiting time distribution for such events.

Fig. 5

Waiting time distribution for local event (f \ge AL^{1/2}) for varying A and fixed (L=256, , 400.)

In Fig. 5a and 5b, the y axis shows the normalised waiting time distribution w(t; A, L) for the interval t between successive large events. This is scaled by a factor of (L^{1/2}), and is plotted against the scaled waiting time (tL^{-1/2}). Notice that the graphs all begin from (0, 0) and peak at a time (t_m). In each graph, the peak (mode), (t_m) is independent of A. The height of this peak falls with A, because the distributions become increasingly long tailed for large A.

3.2.2 Average time spacing between events

In Fig. 6 we plot the log of scaled average waiting time, (\ln \left( \overline{t}L^{-1/2} \right) ) against A, using the Monte Carlo data for (L=256,400). There is a simple argument to get the dependence of (\overline{t}) on A.

From the definition of the average time spacing between such ‘large’ events, we have (\overline{t}^{{-1} = {\text {Prob}},(f \ge AL}{1/2})). Avalanches in the lattice will result in outflow throughout the bottom layer (x \in {1,..,mL^{1/2}}) , but we are interested in the outflow in a smaller interval (1 \le x \le L^{1/2}). Let the total outflow in the bottom layer be (f_g), and that within our window be f. Then, in the large L limit, the probability density for the scaled drop number (v :=fL^{-1/2}) is labelled as h(v). With this,

$$\begin{aligned} \overline{t}(A)^{{-1} = L}{-1/2}\int _A^{1}dv , h(v) , . \end{aligned}$$

(2)

It is a well known result for the 2-neighbour directed sandpile that the avalanche cluster has no holes. So, the two boundaries of the avalanche constitute an annihilating random walk. Suppose a particle is added at the top of the lattice at a co-ordinate (r_0). Consider the annihilating random walk on the lattice corresponding to the ensuing avalanche. Let (x_1 < x_2) be the co-ordinates of the two annihilating particles, and define (x = x_2-x_1 , , r = (x_1+x_2)/2). In these co-ordinates, the conditional density (\rho (x,r;y,|, 1,r_0;0)) evolves as:

$$\begin{aligned} \partial _y\rho (x,r ;y) = \frac{1}{2}\left[ \partial _r^2 + 4\partial _x^2 \right] \rho (x,r;y) , . \end{aligned}$$

(3)

Further, if we integrate out the centre of mass co-ordinate r in (3), we get the probability density p(x; y|1, 0) for a Brownian excursion [24]. The quantity x(L) is the number of toppled sites in the bottom layer, so (x(L) = f_g). A well known result for Brownian excursions is that (p(x,y=L|1,0) = L^{-1/2}g_{b}\left( xL^{-1/2}\right) ), where (g_b(z) = \left( 4\pi \right) ^{{-1/2}z e}{-z^2/4}). Define the intervals on (\mathbb {R}_{mL^{1/2}}) : (I_0 = (0,L^{1/2}]) and (I_{x,r} = \left[ r-\frac{x}{2} , r + \frac{x}{2} \right] ). Then we have:

$$\begin{aligned} h(u) = \int _{0}^{mL{1/2}} \frac{dr_0}{mL^{1/2}} \int _0^{mL^{1/2}} dr , dx , \rho (x,r;L|1,r_0,0) \times \delta \left( uL^{1/2} - |I_{r,x} \cap L^{1/2}I_0| \right) , . \end{aligned}$$

(4)

The first factor of (1/mL^{1/2}) is because (r_0) is uniformly distributed. Using translation invariance, (\rho (x,r;L|1,r_0,0) = \rho (x,r-r_0,L|1,0,0)). So, we get (after also rescaling all quantities by (L^{1/2}), and defining the scaled interval (\tilde{I}_{r’,x’} = \left[ r’ - \frac{x’}{2} , r’ + \frac{x’}{2} \right] ) on (\mathbb {R}_m)):

$$\begin{aligned} h(u) = \frac{1}{m} \int _0^m dx’ , g_b(x’) \int _0^m dr’ , \delta \left( u - |\tilde{I}_{r’,x’} \cap I_0| \right) , . \end{aligned}$$

(5)

These integrals are easily evaluated. We get:

$$\begin{aligned} h(u) = \frac{e^{-u^2/4}}{2\sqrt{\pi }m}\left[ 4 + u(1-u) \right] + \frac{\delta (u-1)}{2 m} ,{\text {erfc}}\left( \frac{1}{2}\right) . \end{aligned}$$

(6)

Finally, using (2),

$$\begin{aligned} \overline{t}(A) = mL^{1/2}\left[ {\text {erfc}}\left( \frac{A}{2} \right) -\frac{1}{2}{\text {erfc}}\left( \frac{1}{2}\right) +\frac{A-1}{\sqrt{\pi }}e^{-A^2/4}\right] ^{-1} , . \end{aligned}$$

(7)

Fig. 6

Average waiting time vs A, for (L = 256, , 400)

3.2.3 SROT curves for outflow variables

Here we discuss MC results and plot SROT curves. The plots are generated using the time series for (L = 144,256,400) and a multiple values of A (Fig. 7a, 7b).

Fig. 7

SROT plots. (a) Plots for (0.5 \le A \le 0.8) (b) Curves are identical for (L = 144, , 256, , 400) with (A = 0.5, , 0.8)

We make the following observations. Firstly, for a given A, there is a scaling collapse in the SROT plots for a very wide range of L. At large lengths, there are only small deviations from the true thermodynamic limit; therefore there is finite predictability in the thermodynamic limit.

Next, there is an approximate ‘on time’ fraction (\nu ^*) such that for (\nu > \nu ^{*), c is close to 1. The value of (\nu *) increases with A. Plots for larger A are similar to those for smaller A, except that they are ‘pushed’ to the right. To understand this we examine the waiting time distribution (Fig. 5b). Fix L and A. The function w(t, A, L) is peaked about (t_m = t_m’L{1/2}). If the alarm is switched off for a time (\tau < t_m), then the average number of large events occurring in the of interval is pretty small, and nearly all avalanches are predicted. This corresponds to (\nu > \nu _m}*) in Fig. 7b. Once we cross the peak, i.e. (\tau >t_m) and (\nu < \nu ^*), there is a surge in the number of missed avalanches and the predicted fraction suddenly drops. On increasing A, the height of the peak, (w(t_m)) decreases and w becomes more and more heavy tailed. Suppose we keep the alarm off for a fraction (\tau < t_m) of the time. Then, the alarm is switched on until the next large event. Due to the heavy tailed nature of w, the next large event takes a much longer time to occur on average. So the on time fraction is much larger, and the ROC curve is pushed to the right for larger A. This means that very extreme events ((A \sim 1)) are in this sense ‘less’ predictable, which is different from other results [14, 15], but not contradictory. The result is specific to the current strategy. In fact, for subsequent strategies, we will see that predictability improves with rarer and larger events. Moreover, it is also true that for (\nu > \nu ^*), the c is larger for larger A, though all c values are close to 1.

A possible extension of this strategy involves using the conditional waiting time distributions given the size of the previous large event. This may improve the prediction efficiency.

3.2.4 What are the SROT curves ?

Under the assumption that the time intervals between successive large events are i.i.d variables, and that their distribution w(t) is known, we can calculate the SROT curves. Fix the alarm-off threshold (\tau ). Suppose that an event happens at (t_0 = 0). There can be multiple ‘missed’ events in the interval ((0,\tau )). Suppose there are n of them, at times (t_1\le t_2 \le ..\le t_n). Let (t_{n+1} > \tau ) be the next ‘detected’ event. Since the time evolution is statistically stationary, it suffices to consider the statistics between two detected events. The c and (\nu ) parameters (true positive fraction and alarm on time fraction) for the SROT curves can be written as :

$$\begin{aligned} c(\tau )&= \sum _{n=0}^{{\infty } \frac{1}{n+1} \int _{\tau }}{\infty } dt_{n+1} \times {\text {Prob}}(n {\text { events in }} (0, \tau ) , | , n+1 ^{{\text {th}} {\text { event at }}t_{n+1} ) \nonumber \ \nu (\tau )&= 1 - \tau \sum _{n=0}}{\infty } \int _{\tau }^{\infty } \frac{dt_{n+1}}{t_{n+1}} \times {\text {Prob}}(n {\text { events in }} (0, \tau ) , | , n+1 ^{\text {th}} {\text { event at }}t_{n+1} ) , . \end{aligned}$$

(8)

Let (\tilde{c}(r) := \int _{0}^{{\infty } d\tau e}{-r\tau }c(\tau )) and (\tilde{\nu }(r) := \int _{0}^{{\infty } d\tau e}{-r\tau }\nu (\tau )) be the Laplace transforms of (c(\tau ), \nu (\tau )). From the i.i.d property we have :

$$\begin{aligned}&{\text {Prob}}(n {\text { events in }} (0, \tau ) , | , n+1 ^{{\text {th}} {\text { event at }}t_{n+1} ) \nonumber \&\quad = \int _{0}}{\tau } dt_1 \int _{t_1}^{{\tau } dt_2 , ..., \int _{t_{n-1}}}{\tau } dt_{n} , w(t_{n+1}-t_n)\prod _{j = 0}^{n-1} w(t_{j+1} -t_{j}) , . \end{aligned}$$

(9)

Let (q(t_n) := {\text {prob}} (n-1 {\text { events in }} (0,t_n) )) where the (n^{\text {th}}) event happens at (t_n). Then,

$$\begin{aligned} q(t)&= \int _{0}^{{t_n} dt_1 .., \int _{t_{n-2}}}{t_n} dt_{n-1} , \prod _{j = 0}^{n-1} w(t_{j+1} -t_{j}) \Rightarrow \tilde{q}(s) = \left[ \tilde{w} (s)\right] ^n , . \end{aligned}$$

(10)

Above, (\tilde{q}(s)) is the Laplace transform of q(t). Let us first evaluate c.

$$\begin{aligned} \tilde{c} (r)&= \int _0^\infty d\tau , e^{-r\tau } \sum _{n=0}^{{\infty }\frac{1}{n+1} \int _{\tau }}{\infty } dt_{n+1} \int _{0}^{\tau } dt_n , q(t_n) , w(t_{n+1} - t_n) \nonumber \&= \frac{1-\tilde{w}(r)}{r} \times \int ^* \frac{ds}{2\pi i} \frac{\ln \left[ 1-\tilde{w}(s) \right] }{\tilde{w} (s)(r-s)} , . \end{aligned}$$

(11)

In the second line we have used a Mellin integral for the inverse Laplace transform. Next, we evaluate (\tilde{\nu }(r)).

$$\begin{aligned} \tilde{\nu }(r)&= \int _0^{\infty }d\tau , \tau e^{-r\tau }\sum _{n=0}^{{\infty }\int _\tau {\infty } \frac{dt_{n+1}}{t_{n+1}} \int _0{\tau } dt_n , w(t_{n+1}-t_n) , q(t_n) \nonumber \&= \frac{1}{r}\int * \frac{ds}{2\pi i} \frac{1}{1-\tilde{w}(s)} , \left[ \frac{\tilde{w} (r)}{s-r}-\int _{r}{\infty } \frac{dk , \tilde{w}(k)}{r(k-s)} + \int _0}\infty \frac{dk , \tilde{w}(k), (k+2r-s)}{r(k+r-s)^2}\right] , . \end{aligned}$$

(12)

3.2.5 Special limits

Next we check the limiting behaviour of the SROT plots, i.e. the top right and bottom left regions of Fig. 7a. These correspond to (\tau \ll t_m) and (\tau \gg t_m) respectively, where (t_m = t_m’L^{1/2}) is the mode of w(t). Let (\tau = \tau ’L^{1/2}).

In the first case, the alarm is switched off only for a short interval corresponding to the (t’ \ll t_m’) region in Fig. 5a. Let (w\left( t’L^{1/2}\right) _{t’ \ll t_m’} \sim w_0(t’)). For instance, it may be a power law. In the short off interval, the probability of n missed events falls very quickly with n, so we can keep only the first few small n terms in the calculation of c and (\nu ).

$$\begin{aligned} c(\tau ) \simeq \int _{\tau }^{{\infty } , dt_1 w(t_1) +\frac{1}{2} \int _{\tau }}{\infty } dt_2 \int _{0}^{{\tau } dt_1 , w(t_2-t_1)w(t_1)&\simeq 1-\frac{1}{2}\tau w_o(\tau ) . \nonumber \ 1-\nu (\tau ) \simeq \tau \int _{\tau }}{\infty } dt_1 , t_{1}^{-1} w(t_1) + {\text {H.O.}}&\simeq \ \frac{\tau }{t_m} = \frac{\tau ‘}{t_m’}, . \end{aligned}$$

(13)

With (\overline{c} = 1-c , , \overline{\nu } = 1-\nu ), the small (\tau ) behaviour of the plot is :

$$\begin{aligned} \overline{c} \simeq \frac{1}{2} \overline{\nu }, t_m’ , L^{-1/2}, w_0\left( \overline{\nu }, t_m’\right) , . \end{aligned}$$

(14)

For instance, if (w_0(t’) \sim t’^{{\gamma }), (\overline{c} \sim \overline{\nu }}{\gamma +1}). Note that in the thermodynamic limit, (\overline{c} \rightarrow 0).

In the opposite limit, the alarm is switched off for a very large time (\tau \gg t_m). The distribution for the number of avalanches in the interval ((0, \tau )) is sharply peaked about (n^* = \tau /t_m). (This can be shown by taking a Gaussian approximation of w(t) about (t_m), in which case only the events around the peak (t_m) dominate the probability measure). We have :

$$\begin{aligned} c&\simeq \frac{t_m}{\tau } , , , \nu \simeq 1-\tau \int _0^{\infty } \frac{ dt , w(t)}{t + \tau } \simeq \frac{\overline{t}}{\tau } , . \nonumber \ c&\sim \frac{t_m}{\overline{t}} \nu , . \end{aligned}$$

(15)

Above, (\overline{t} \equiv \overline{t}(A)) is the average spacing time given in (7). For increasing A, we know that (\overline{t}(A)) increases, whereas (t_m) is independent of A. One can see in the plots (Fig. 7b) that for larger A, the slope of the curve near (0, 0) is diminished.

3.2.6 Other local variables : toppling numbers and heights.

The previous sections have dealt with drop number variables. But the equations (11), (12), (14) and (15) were derived for an arbitrary waiting time distribution w(t). So, they hold good also for toppling number and height events. Of course, we must use the appropriate waiting time distributions (w^{S}(t)) and (w^h(t)) in the two cases. From the time series data, we plot the local waiting time distributions for toppling number events (S \ge AL^{3/2}), with (L=256) and a range of A (Fig. 8). Similar to the outflow curves, the peak in these curves is also independent of A.

Fig. 8

Local waiting time distribution for (S \ge AL^{3/2})

4 Strategies using the average excess density in the pile

The sandpile model involves a Markov update rule for the full configuration. The time series ({ S_t} , {h_t} {\text { and }} {f_t}) constitute Hidden Markov processes [25]. To enhance the predictability of the observables, one has to incorporate some more information about the sandpile configuration. We will use the excess density (\epsilon \in (-1,1)), which is defined by the total mass on the lattice being (N(1+\epsilon )/2). We are given a past-time series (\epsilon _1, \epsilon _2,..,\epsilon _t). We want to predict events of the form (S \ge AL^{\gamma }), (h \ge zL^{\beta }) or (f \ge AL^{\delta }) where (z, A \sim \mathcal O(1)), and (\gamma , \beta ,\delta ) are appropriate scaling exponents. The strategy fixes the prediction variable as :

$$\begin{aligned} X_{t+1} = {\left{ \begin{array}{ll} 1 , & \epsilon _t \ge \epsilon ^* \ 0, & \epsilon _t < \epsilon ^* , . \end{array}\right. } \end{aligned}$$

(16)

Two trivial cases can be immediately understood. In the case (\epsilon ^* \rightarrow -1), the alarm is always on. Consequently, all avalanches are predicted and (c = \nu = 1). Conversely, for (\epsilon ^* \rightarrow +1) the threshold is too high and the alarm is almost always off. So (c = \nu = 0). We can vary the parameter (\epsilon ^*) in ([-1,1]) to get SROT curves (c(\epsilon ^*)) vs (\nu (\epsilon ^*)).

4.1 Height variables

We discuss the behaviour of predictions for events of the form (h \ge zh^*).

4.1.1 Scaling theory

We have the unscaled variables (h^*, \epsilon ). Let us define scaling exponents so that scaled variables are (z := h^*L^{-\alpha }, , \xi := \epsilon L^{\beta }). Further, set (M = mL^{\mu }). In this strategy, the prediction for time (t+1) is blind to all data other than that at t. So, one must average over all possible histories of the system and the measure that we use is same as the steady state measure. The x and y axis of our SROT curve are given by :

$$\begin{aligned} c&= \frac{{\text {Prob}} \left( h \ge zL^{\alpha } , ; , \epsilon \ge \xi L^{-\beta } \right) }{{\text {Prob}} \left( h !\ge ! zL^{\alpha }\right) } !\equiv ! \frac{1}{{\text {Prob}} \left( h !\ge ! zL^{\alpha }\right) }! \int _{\xi L^{-\beta }} {\text {Prob}}(d \epsilon ’) !\times ! {\text {Prob}} \left( h! \ge !zL^{\alpha } ,| , \epsilon ’ \right) , . \end{aligned}$$

(17)

$$\begin{aligned} \nu&= {\text {Prob}} \left( \epsilon \ge \xi L^{-\beta } , \right) , . \end{aligned}$$

(18)

All integrals are over the steady state measure ({\text {Prob}}(d \epsilon )). In the steady state all states are equally probable. This is due to the fact that in the directed model, all states are recurrent [26]. So, the probability to have an excess mass (N\epsilon /2) in the system is simply (2^{-2N}\left( {\begin{array}{c}N\ \frac{N}{2} (1+\epsilon )\end{array}}\right) ). Here, (N := L^{1+\mu }) is the total number of lattice sites. In the thermodynamic limit, the density is:

$$\begin{aligned} \rho \left( \epsilon = \xi L^{-\beta }\right) \sim \frac{m^{-\frac{1}{2}}L^{-\frac{1+\mu }{2}}}{\sqrt{2\pi }}\exp \left( -\frac{m\xi ^{2}{2} L}{1+\mu -2\beta }\right) , . \end{aligned}$$

(19)

The typical density fluctuation is thus of the order \