Sandpile-based model for capturing magnitude distributions and spatiotemporal clustering and separation in regional earthquakes

We propose a cellular automata model for earthquake occurrences patterned after the sandpile model of selforganized criticality (SOC). By incorporating a single parameter describing the probability to target the most susceptible site, the model successfully reproduces the statistical signatures of seismicity. The energy distributions closely follow power-law probability density functions (PDFs) with a scaling exponent of around − 1.6, consistent with the expectations of the Gutenberg–Richter (GR) law, for a wide range of the targeted triggering probability values. Additionally, for targeted triggering probabilities within the range 0.004–0.007, we observe spatiotemporal distributions that show bimodal behavior, which is not observed previously for the original sandpile. For this critical range of values for the probability, model statistics show remarkable comparison with long-period empirical data from earthquakes from different seismogenic regions. The proposed model has key advantages, the foremost of which is the fact that it simultaneously captures the energy, space, and time statistics of earthquakes by just introducing a single parameter, while introducing minimal parameters in the simple rules of the sandpile. We believe that the critical targeting probability parameterizes the memory that is inherently present in earthquakegenerating regions.


Introduction
The sandpile model, introduced as a representative system for illustrating self-organized criticality (SOC; Bak et al., 1987), has opened up new avenues for the use of discrete cellular automata (CA) models in capturing the salient fea-tures of many systems in nature (Olami et al., 1992;Drossel and Schwabl, 1992;Malamud and Turcotte, 2000;Piegari et al., 2006;Juanico et al., 2008).Seismicity, which is rife with power-law statistical distributions (Saichev and Sornette, 2006), is an interesting test case for such approaches.Despite the complexity of the processes in the earth's crust that limit our ability for accurate, short-term prediction of events, it is worth noting that many statistical features of seismicity, as obtained from substantially complete earthquake records, can be recovered using simple CA models.
One of the earliest attempts for sandpile-based modeling of earthquake distributions is by Bak and Tang (1989), who used a two-dimensional sandpile to show the powerlaw Gutenberg-Richter (GR) distributions of earthquake energies (Gutenberg and Richter, 1954).Subsequent authors also noted that the simple sandpile produces power-law distributions of earthquake waiting times upon introducing a threshold magnitude (Paczuski et al., 2005).Additional parameters have been introduced in the model to account for other features of seismicity.Ito and Matzusaki (1990) introduced aftershock triggering to the sandpile model to recover the aftershock frequencies and the hypocenter distributions, which also follow power-law decays.To represent a scale-invariant distribution of earthquake faults, Barriere and Turcotte (1991) incorporated a power-law distribution of box sizes in the CA model and recovered not only the GR distribution but the occurrence of foreshocks.On the other hand, Steacy et al. (1996) investigated the effect of a heterogeneous strength distribution and found that the power-law exponent of the magnitude distribution is dependent on the degree of the heterogeneity.Inspired by the sandpile design, Olami et al. (1992) used a CA implementation of the earlier Published by Copernicus Publications on behalf of the European Geosciences Union & the American Geophysical Union.
Burridge-Knopoff model (Burridge and Knopoff, 1967) that incorporates dissipative terms and inhomogeneous energy redistribution rules to capture key elements of seismicity, along with foreshocks and aftershocks (Hergarten and Neugebauer, 2002).In another work, Jagla (2013) has shown that the GR law can be recovered in a forest-fire model, with the fires interpreted as the earthquake occurrences.
The introduction of additional parameters to subsequent models indicates that the simplest rules of the original sandpile are not able to capture key features of seismicity.In the sandpile model, the stress in the grid is released in a single avalanche event resulting from small-neighborhood cascades; for seismicity, the energy is released in a sequence of correlated events.Additionally, the single triggering at random locations will tend to produce normal distributions of interoccurrence distances and times, which, again, deviate from those observed in records of seismicity.Finally, the conservative sandpile with symmetric nearest-neighbor redistribution rules does not take into account the memory that may be present in actual earthquake-generating zones.
In this work, we adhere to the key features of the sandpile model, and introduce a very simple modification: for a fraction of the iteration times, determined randomly, we direct the triggering into the most susceptible site in the grid.In this case, the avalanches in the grid are deemed to be analogous to the energy release during an earthquake occurrence.Interestingly, this very simple modification in the sandpile rule enabled us to recover, simultaneously, the distributions of event sizes, interevent distances, and interevent times that are comparable to those obtained from substantially complete earthquake records.

Model specifications
The model utilizes a two-dimensional space discretized into a grid of L × L cells arranged in a square lattice.The cells contain continuous-valued information states σ representing the local measure of susceptibility to rupture.At time t = 0, the states are initialized to have values within [0, σ max ), where, in this case, we set σ max = 1.0 as the relative measure of the rupture threshold.
The dynamical evolution of the grid is guided by rules patterned after the Zhang sandpile that uses continuous-valued states (Zhang, 1989).We choose an asynchronous update rule, such that every discrete time step the grid is triggered by adding a constant value ν to a single location (x, y), σ (x, y, t + 1) → σ (x, y, t) + ν, with all the other sites unperturbed.The asynchronicity may represent the nonuniformity of crustal motion that drives the accumulation of elastic potential energy at faults.Moreover, the model introduces a targeted triggering probability p that the most susceptible site, i.e., the site with the highest σ value in the grid, will receive the driving term ν.Triggering is therefore applied to the most susceptible site with probability p and to a randomly chosen site with probability 1−p.The value of p represents a memory term, and parameterizes the tendency of fracture to occur at more susceptible locations along an earthquake-generating zone.
In the event that a cell matches or exceeds a maximum possible value σ max , the local region is deemed to rupture.No new trigger is added to the system during such events; instead, the stress from the collapsing site is transferred to the four nearest neighbors in the grid, σ (x ± 1, y ± 1, t) → σ (x±1, y±1, t)+ 1 4 σ (x, y, t), leading to the relaxation of the original site, σ (x, y, t) → 0. Such relaxations may produce a cascade of subsequent stress redistributions and relaxations in the grid when one or more of the neighbors are driven to the threshold.As in the previous sandpile models, the number of affected sites in the grid, A, is tracked to quantify the relative event size.Additionally, we also recorded the number of unique activations V , the number of times a cell has been affected by a cascade, as a proxy for the actual energy or seismic moment of the relaxation event.
Prior calibrations show that ν = 10 −3 produces power-law event size distributions comparable to the GR law, and that t max = {1, 4, 16} × 10 7 iterations, where the first 10 % are neglected for transient behavior, produces a substantial number of avalanche events for L = {256, 512, 1024} grids, respectively.We investigated the case of different targeted triggering probabilities, p = {0, 1 × 10 k , 5 × 10 k , 1}, where the integer k ranges from −5 to −1, to scan a wide range of possible system behaviors.For each of the p values, we track all nonzero A i and V i and their avalanche origins and occurrence times (x i , y i , t i ), where i denotes the temporal index of occurrence of an event.The spatial and temporal separations of successive events, R i = [(x i −x i−1 ) 2 +(y i −y i−1 ) 2 ] 1/2 and T i = t i −t i−1 , are computed, and the probability density functions (PDFs) of all A, V , R, and T values are plotted.
Records of very low-magnitude earthquakes are oftentimes incomplete because they are both too weak for detection and their occurrence is orders of magnitude in frequency as compared with the higher-magnitude ones.In the model, however, we can resolve all the avalanches, even the smallest ones that affect only single neighborhoods.To mimic the effect of the non-retention of the smallest earthquakes, we employed a thresholding procedure in the analyses by setting A th = {5, 10, 50, 100, 500, 1000, 5000} such that all events with A < A th are removed from the sequence.Because the A PDF is just expected to be cut off below A th , we observed how the statistical distributions of R and T will be affected upon employing different A th values.
Finally, as a way of comparison and verification, we compare the model statistics with those obtained from actual earthquake catalogs from Japan (JP), Philippines (PH), and southern California (SC), as investigated in a previous work by Batac and Kantz (2014).The JP records are obtained from the Japan University Network Earthquake Catalog (JUNEC), with approximately 137 000 events from July 1985-December 1998; the PH earthquakes are composed 10 1 10 2 10 3 10 4 10 5 A (cells) 10 -6 10 -3 10 0 Prob(A) 10 -14 10 -8 Prob(E) of 70 000+ events from 1973 to 2012, as obtained from the Preliminary Determination of Earthquakes (PDE) Catalog; while the SC records are from the Southern California Earthquake Catalog (SCEC) containing 516 000+ events from 1982 to 2012 (events due to man-made activities are removed).We compared the behaviors of the model and data statistics using scaling factors derived from model parameters.

Model results
Figure 1a shows the avalanche size PDFs for the different values of the targeted triggering probability p.For the broad range of p values considered, the distributions are found to be comparable to a power-law A −α with α = 1.6.Continuous-state sandpiles have been known to have avalanche size scaling exponents greater than 1.0, the exponent of the discrete Bak-Tang-Wiesenfeld (BTW) sandpile.Lübeck (1997) conducted large-scale simulations of a similar Zhang sandpile and obtained exponents slightly higher than 1.2, which can go even higher for large driving rates ν.In a similar model that incorporated non-conservation, Piegari et al. (2006) obtained power-law exponents that approach 1.6 in the conservative limit for the same order or magnitude of ν that we used.The higher exponents and the effect of the driving rates are also verified by an equivalent conservative model and actual sand avalanche experiments by Juanico et al. (2008), and in other asynchronous updating models (Paguirigan et al., 2015).The resulting power-law exponent is deemed to be a result of the accumulation of stress at various locations; be-cause the triggering is done at only a single site every time, there is little global connectivity among critical sites, resulting in a preponderance of smaller, isolated avalanches.The fact that the distributions are almost similar regardless of the value of p indicates that the targeted triggering probability has minimal effect on the avalanching mechanism of the grid, such that the system preserves the SOC characteristics of the original sandpile.In contrast, the Olami-Feder-Christensen (OFC) model, one of the foremost discrete models of seismicity, tends to lose the universality of the exponents upon the introduction of non-conservation (Olami et al., 1992).
In Fig. 2a, we observe that the original sandpile p = 0 produces unimodal statistics, whose tails decay towards the largest possible distance √ 2L in the finite grid.The simple sandpile, therefore, is not capable of replicating the observed earthquake separation distance distributions, which are found to exhibit bimodality due to the difference in the characteristic times of the correlated aftershock sequences and the independent mainshocks (Baiesi and Paczuski, 2004;Zaliapin et al., 2008;Touati et al., 2009;Batac and Kantz, 2014).This inspired the introduction of p, which is a random occurrence in time but is inherently affecting the spatial distribution of events in the grid.We do note here that the parameter p is just the probability to target the most susceptible site in the lattice, unlike previous implementations that actually preselect the next targeting location within the vicinity of the previous avalanche (Ito and Matzusaki, 1990).Indeed, without the imposition of such a spatial bias, the replication of the short-R regimes is not guaranteed.Interestingly, however, the plots in Fig. 2a show increased probability of occurrence of the short-R distances upon introducing nonzero p. From this, we can deduce that the most susceptible sites in the lattice are most likely to be found within the vicinity of a previous large avalanche, a fact that was not exploited by earlier similar models.In fact, in the biased case p = 1, we recovered unimodal statistics, as shown in Fig. 2a, albeit at a shorter characteristic distance; for the L = 256 grid, the average location of the most susceptible site from the previous avalanche origin was found to be around 21 cell lengths.
Midway between these two extremes (p = 0 for the original and p = 1 for the completely biased sandpile), we can find a suitable value of p where reasonable comparison with empirical data can be obtained.
The interevent time distributions are shown in Fig. 3a for the L = 256 and t max = 10 7 iterations.We observe the expected shift of the tail cutoff towards shorter T values as p is increased; triggering the highly susceptible sites will more likely result in a new avalanche event, thereby shortening the average waiting time.The resulting distributions are for the case wherein all the events are included in the sequence; we expect a lengthening of the tails of the distributions when we neglect other events below the threshold A th .

Energy distributions and the Gutenberg-Richter law
The GR law, which is usually presented in terms of the magnitude m and as a complementary cumulative distribution function (CCDF) log 10 m = a − bm, can be shown to be equivalent to an energy E CCDF that behaves as E −2/3 from the definition of m and by assuming b = 1, which is the case for most complete records (Jagla, 2013).By noting that the CCDF is effectively an integral of the PDF, the earthquake energy PDF will then behave as E −5/3 .In Fig. 1b-d, similar power-law trends have been obtained for the JP, PH, and SC records, which have different levels of catalog completeness, as indicated by the extent of the power-law regimes.To minimize the problems associated with the inherent incompleteness of smaller-energy events (Zaliapin and Ben-Zion, 2015), we impose a threshold magnitude m th for succeeding analyses such that earthquake events with magnitudes lower than m th are dropped from consideration.The range of such magnitudes considered, which are well within the power-law regimes of the plots, is shaded in Fig. 1b-d: m th ∈ [2.5, 3.5] for JP and SC and m th ∈ [4.5, 4.8] for PH.
In keeping with the earlier sandpile-based approaches where the avalanche size A is used for comparison with earthquake energies (Bak and Tang, 1989;Ito and Matzusaki, 1990), we present in Fig. 1 the PDFs of A with those of E from the seismogenic regions considered.It is worth emphasizing that similar power-law trends result from the introduction of the parameter p, regardless of how large its relative value is.We note, however, that aside from the avalanche size A, there are other parameters that can be used to track the extent of the avalanche event.One such measure is the number of activations V , wherein the sites repeatedly affected by the avalanching process get to be counted multiple times.Previous works have shown that V and A in discrete models may in fact have actual associations with the seismic moment and fracture area, respectively, and may exhibit nontrivial scaling relations (Landes and Lippiello, 2016).We present in Fig. 4a the distributions obtained upon tracking V .triggered sandpile) results in a V (A) ∝ A 1.5 scaling.On the other hand, for p = 1 (sandpile with targeted triggering), the behavior appears to shift towards V (A) ∝ A 1.3 for very large A values.This lower scaling exponent of the activation for large avalanche sizes is expected for targeted triggering; because the most susceptible site is always targeted, there is minimal accumulation of near-critical sites near the location of the avalanche origin, which results in a lower number of reactivations of affected sites near an avalanche event.

Spatial separation of earthquake events
In the original asynchronous sandpile models, one only recovers unimodal statistics for interevent distances.This is due to the stochastic nature of the triggering: the next location to be perturbed is drawn from an oftentimes uniform distribution; i.e., all sites are likely to be triggered next.Additionally, the nature of internal cascading within the sandpile grid results in the depletion of all the critical sites within the extent of the avalanche area.The same cannot be said of earthquakes: after the release of elastic potential energy at a fault location, the subsequent crustal motion may tend to favor other fractures near the vicinity of the earlier event to release the remaining stored energy.
Interestingly, the addition of the simple targeted triggering probability p has enabled us to recover statistical distributions that are comparable to those observed in regional earthquake records up to a scaling factor.It should be noted that without any form of spatial clustering, the characteristic separation distance is limited by the finite system size.Rescaling is therefore conducted by comparing the characteristic sizes (modes) of the memoryless cases of the model (p = 0) and the data (shuffled sequence).The interevent distance distributions of the shuffled sequences are shown as the black symbols in Fig. 5, while the corresponding model p = 0 distribution is shown in Fig. 2a, with both clearly showing unimodal statistics.
Upon getting the rescaling factor, we scan through the possible p values to obtain p values that will result in comparable R distributions between model and data.We observe that the model parameters that will correspond to the empirical distributions upon such a simple rescaling range from p * ≈ 0.004 to 0.007.Figure 2b-d show the interevent distances between successive earthquakes in the different regional records considered, superimposed with the rescaled statistics of the model.
The rescaled model statistics for p = 0.007 show good agreement with interevent distances from the three seismogenic regions.As expected, larger grid sizes will result in a better discrimination of shorter R; i.e., one pixel unit will correspond to shorter actual distance units.In our case, for the largest grid size used (L = 1024), we find that the scaling factors obtained by matching the modes result in the following correspondence with a unit cell length: 1.3 for JP, 1.2 for PH, and 0.5 km for SC.The distributions are found to be similar regardless of the threshold magnitude A th considered due to the finite system size; even upon removing the weakest events, the avalanche origins are confined within the grid, resulting in the same probability density distribution of R.

Temporal separation of earthquake events
The temporal separation of aftershocks and mainshocks that have different characteristic waiting times is an intuitive result that is both well known and widely studied (Zaliapin et al., 2008;Touati et al., 2009;Batac and Kantz, 2014;Batac, 2016).The proposed model, therefore, must also show these features to be able to compare reasonably well with the temporal distributions of seismicity.In the following, we compare the results of the model having p * = 0.007 and grid dimension L = 1024, which has been shown to have comparable R statistics with empirical data.
In comparing model and empirical temporal interevent statistics, one does not have the similar advantage of having a finite "space."The goal of rescaling in time is to re- cover the relatively short T regimes first; theoretically, the longest T will be recovered if the model is allowed to run for very long iteration times.Additionally, in rescaling the time, one should take into account the fact that the earthquake record is thresholded by m th , effectively lengthening the average time between the occurrence of two events.Ideally, if all the events, no matter how weak, can be detected and recorded, we would not have long tails in the waiting time distribution of earthquakes.This is also observed in sandpile-based models; previous approaches have shown that the waiting time distribution will be Poisson distributed when all the events are considered, but will begin to show apparent power-law characteristics upon thresholding (Paczuski et al., 2005;Juanico et al., 2008).
For our purpose, we arbitrarily chose the following threshold avalanche sizes for removing weaker events: for comparison with JP and SC, which are both taken to have m th = 2.5, we used A th = 5 × 10 3 ; on the other hand, for PH, with relative completeness beyond m th = 4.5, A th = 5 × 10 5 is used.The values of A th are obtained by maintaining the fraction of events left after neglecting the weaker events.Still, because of the limited number of regional data sets considered, which does not allow for further testing their correspondence, we emphasize that the values of the A th obtained do not necessarily translate into an exact equivalence with the threshold magnitude m th for the data.
Upon removing the events with A < A th , we obtained the modes of both the data and the model for visual comparison.This resulted in slight differences in the rescaling factors for the different data sets.One iteration of the model corresponds to 0.006 s for JP, 0.004 s for PH, and 0.002 s for SC.Apart from recovering the qualitative trends in the T PDFs, we conducted additional analyses to check if the model results also show spatiotemporal clustering and separation behaviors.In Fig. 5, we mark the location of R * , the characteristic separation distance where the empirical distributions and those of the shuffled sequences begin to show comparable trends.The R * values of the region considered, which are similar to the results of Batac and Kantz (2014), are deemed to be a good marker for separating "nearby" and "far away" events.We note that a similar procedure done using the rescaled model statistics results in compa- As shown in Fig. 6a-c, for all the seismogenic regions considered, the distributions of T in and T out differ significantly from that of the total T .The relative frequency plots of T in all cases can be shown to be a crossover between T in and T out that have different modes.As expected, the T out distributions do not coincide due to the different periods involved in the catalogs considered.On the other hand, the T in distributions all have modes at short T values, suggesting a strong dependence among the interevent properties in space and time (Livina et al., 2005).This conditional distribution therefore quantifies the spatiotemporal clustering observed in earthquakes, particularly among aftershock sequences that result from the correlated mechanisms: "nearby" events are also more likely to be separated by shorter waiting times.
In Fig. 6d-f, we observe that despite the shorter iteration times being considered, the model was able to show the separation of the T in and T out distributions, a feature that is also found in empirical data (Batac and Kantz, 2014) and in other earthquake models (Touati et al., 2009).Moreover, it is particularly interesting to note that the rescaled T in statistics of model and corresponding T in from the earthquake data show comparable trends, especially for shorter waiting times, as shown in the insets.The T in statistics have been shown to correspond with the statistics of aftershocks, as shown in studies of fresh aftershock statistics from empirical data (Batac, 2016).This suggests that the correlated mechanisms in actual earthquake systems that produce the T in distributions are also present in the model.

Model advantages and insights on empirical modeling
Introducing the parameter p into the sandpile driving is a straightforward way of incorporating memory into the system.This simple parameter holds a distinct advantage over other models that introduced additional parameters, because it spans a wide range of possible statistical distributions in event size, space, and time, without actually biasing the location of the next triggering event.Being a single parameter, the correspondence between p and actual properties of the earthquake-generating system may be difficult, if not impossible, to ascertain.At best, we may think of p as a combined effect of many different factors on the ground that lead to the preferential triggering of a location.We believe that this parameter, which, for earthquakes, show comparable statistics for the range p * ≈ 0.004-0.007,may be introduced in other sandpile-based models of other events in nature deemed to show self-organized (critical) characteristics.It may be possible to quantify the extent of "memory" of these systems through the value of the parameter p that best replicates their statistical distributions.
Moreover, a deeper analysis of the other regimes of p may lead to a better comparison between the model and other similar protocols.For example, for higher values of p, the model may exhibit extremal dynamics, resulting in more avalanche events due to the tendency to always trigger the most susceptible site.On the other hand, for very low values of p, the dynamics may be comparable to other models that employ uniform loading.Knowing these limits, and establishing how similar and/or different the model is from other discrete models may help put the results in a better context.

Conclusions
In summary, we have presented a simple cellular automata model inspired by the original sandpile model.The model avoids introducing biased rules and instead incorporates a probability of targeting the most susceptible site in the grid, reminiscent of the assumed fracture mechanism of actual earthquake systems.Within a small range of values (p * ≈ 0.004-0.007),we have observed that the model statistics show comparable trends with empirical distributions of earthquake occurrences in energy, space, and time, upon simple rescaling.
The work has also uncovered an important property of the sandpile grid: the most susceptible sites lie within the vicinity of a previous large avalanche event.Previous sandpilebased models that synchronously update all lattice sites, or those that asynchronously update at random locations, are not able to exploit this important property, preventing the possibility of directly modeling earthquakes using the sandpile paradigm.The introduction of such a targeting probability without destroying the sandpile properties may hint at self-organized critical mechanisms at work in the grid.The fact that the simple targeted triggering probability simultaneously recovers these important statistical features of earthquakes is a simple yet novel concept that has not been exploited by previously proposed discrete models based on the sandpile.
Deeper analyses and comparisons with other established models of seismicity may help further establish similarities and differences and put the model results in a better context.Additionally, the parameterization of memory in the form of the targeted triggering probability may be extended to other similar models to possibly capture the statistical distributions of other self-organized (critical) events in nature and society.

Figure 1 .
Figure1.Avalanche size and earthquake energy PDFs.For all figures, lines corresponding to the power-law trend with exponent α = 1.6 are provided as guides.(a) Model results show similar behaviors despite the large differences in p, signifying the retention of sandpile characteristics.The obtained power-law distributions are comparable to the power-law trends in the energy distributions from (b) Japan (JP), (c) Philippines (PH), and (d) southern California (SC).In panels (b)-(d), the horizontal axes scales are preserved; shaded regions denote energy values with substantial completeness which will be used for subsequent analyses.

Figure 2 .
Figure 2. Interevent distance statistics of model, with rescaling for comparison with actual earthquake separation distance data.(a) For an L = 256 grid, higher p results in the preponderance of short-R values.The trends of the model closely mimic those of the data for (b) JP, (c) PH, and (d) SC, where calibration was done by comparing the modes of the model p = 0 and shuffled sequences of the empirical data.Larger grids in panels (b)-(d) result in the capability to replicate the shorter R regimes.

Figure 5 .
Figure 5. Empirical earthquake interevent distance distributions (hollow symbols), along with the corresponding shuffled sequences (filled symbols) for (a) JP, (b) PH, and (SC).The broken lines indicate the R * values where the original and the shuffled sequences begin to show similar trends.

Figure 6 .
Figure 6.Conditional relative frequency distributions of T in and T out for (a)-(c) earthquake data and (d)-(f) corresponding rescaled model results, plotted with the relative frequency plot of all T values.Nearby (far away) events have a higher (lower) chance of having short waiting times and a lower (higher) chance of having long waiting times, as can be seen from the modes of the conditional frequency distributions.The insets of panels (d)-(f) show that the T in PDFs of model and rescaled data have significant overlap, signifying the similarities in their correlated origins.
Fig. 3bd above show the rescaled model distributions alongside the those of the empirical data, showing qualitative similarities in their trends.
rable R * values.Using R * = 164 km for JP, R * = 125 km for PH, and R * = 79 km for SC, we separated the corresponding waiting times T into the sets T in = {T |R ≤ R * } and T out = {T |R > R * }. Figure 6 shows the relative frequency plots of T , superimposed with those of T in and T out , for the empirical data and the rescaled model values.