Recurrent frequency-size distribution of characteristic events

Many complex systems, including sand-pile models, slider-block models, and earthquakes, have been discussed whether they obey the principles of self-organized criticality. Behavior of these systems can be investigated from two different points of view: interoccurrent behavior in a region and recurrent behavior at a given point on a fault or at a given fault. The interoccurrent frequency-size statistics are known to be scale-invariant and obey the power-law Gutenberg-Richter distribution. This paper investigates the recurrent frequency-size behavior of characteristic events at a given point on a fault or at a given fault. For this purpose sequences of creep events at a creeping section of the San Andreas fault are investigated. The applicability of the Brownian passage-time, lognormal, and Weibull distributions to the recurrent frequency-size statistics of slip events is tested and the Weibull distribution is found to be a best-fit distribution. To verify this result the behaviors of the numerical slider-block and sand-pile models are investigated and the applicability of the Weibull distribution is confirmed. Exponents of the best-fit Weibull distributions for the observed creep event sequences and for the slider-block model are found to have close values from 1.6 to 2.2 with the corresponding aperiodicities of the applied distribution from 0.47 to 0.64.


4
. All these systems have been discussed in the literature whether they obey (under some conditions) the principles of SOC.
To estimate hazard risks for a region it is necessary to know the statistical properties of earthquake occurrence. The interoccurrent statistics play an important role in these assessments. Therefore many studies have been devoted to the investigation of the interoccurrent behavior of earthquakes. However, the knowledge of only interoccurrent properties is not sufficient. It is also necessary to know the statistical properties of another type of behavior -the recurrent behavior at a given point on a fault or at a given fault. Unfortunately, about this specific type of behavior we have little information. The reason is that this type of behavior is much more difficult to be investigated. For the interoccurrent frequency-size statistics of a region it is only necessary to count magnitudes of earthquakes occurred in this region. For the recurrent behavior it is required to associate these earthquakes each to a specific fault or fault segment. Actually, the question what statistics correspond to the recurrent behavior of earthquakes has not been answered yet. In the literature primarily the recurrent time-interval statistics have been investigated (Abaimov et al. 2007a;Abaimov et al. 2007b;Hagiwara 1974;Jackson et al. 1995;Matthews et al. 2002;Molchan 1990Molchan , 1991Nishenko and Buland 1987;Ogata 1999;Rikitake 1976Rikitake , 1982Rikitake , 1991Utsu 1984; Working Group on California Earthquake Probabilities 5 1988,1990,2003). Only a few attempts have been made to investigate the recurrent frequency-size statistics (e.g., Abaimov et al. 2007a;Bakun et al. 2005).
This paper overcomes this lack. As some possible alternatives for the recurrent frequency-size behavior we will consider the Brownian passagetime, lognormal and Weibull distributions. In Section 2 we briefly describe these distributions. The tests of goodness-of-fit that we use are described in Section 3. These include the Kolmogorov-Smirnov test and the root mean squared error test.
Ideally, observed sequences of earthquakes on a fault should be used to establish the applicable statistical distribution. However, the numbers of events in observed earthquake recurrent sequences are not sufficient to establish definitively the validity of a particular distribution (Savage 1994).
To illustrate this in Section 4 we will consider the sequence of earthquakes on the Parkfield segment of the San Andreas fault.
In Section 5 of this paper we study sequences with a large number of recurrent events. For this purpose we investigate a creeping section of the San Andreas fault. Creep events on the central section of the San Andreas fault have been studied extensively. The creep measurements have been carried out since the 1960s by U.S. Geological Survey and show both steady-state creep and well defined slip events. We consider the recurrent 6 statistics of slip events that are superimposed on the steady-state creep. For the four sequences we test the applicability of the Brownian passage-time, lognormal, and Weibull distributions to the recurrent frequency-size statistics of slip events. In this case, we have enough events to convincingly differentiate between alternative proposed distributions. In each case we compare the data (sample distribution) with the three distributions and provide tests of goodness-of-fit.
Although the creep records provide enough (up to 100) events in a sequence to differentiate between the proposed distributions, the applicability of these statistics to earthquakes could be questionable.
Therefore, it would be reasonable to test our results with the aid of numerical simulations. A slider-block is often used to study earthquake behavior (see e.g. Abaimov et al. 2007b;Abaimov et al. 2008;Carlson and Langer 1989). Although the applicability of this model to earthquakes is often questionable too, this would give us an independent verification of our results. Therefore in Section 6 we investigate the recurrent behavior of the slider-block model.
A sand-pile model was a basic model for developing SOC ideas (Bak et al. 1988). Therefore it would be curious to see what recurrent behavior this model has. We investigate this question in Section 7. 7

APPLICABLE DISTRIBUTIONS
Three widely used statistical distributions in geophysics are the Brownian passage-time, lognormal, and Weibull. We will compare each of them with our data (sample distribution).

BROWNIAN PASSAGE-TIME DISTRIBUTION
The probability density function (pdf) of amplitudes S for the Brownian passage-time distribution is given by (Chhikara and Folks 1989) where µ is the mean, σ is the standard deviation, and μ σ = V C is the aperiodicity (coefficient of variation) of the distribution. The corresponding cumulative distribution function (cdf) can be obtained analytically (Matthews et al. 2002) but the expression is lengthy and is not given here explicitly.

LOGNORMAL DISTRIBUTION
The lognormal is one of the most widely used statistical distributions in a wide variety of fields. The pdf of amplitudes S for the lognormal distribution is given by (Patel et al. 1976 The lognormal distribution can be transformed into the normal distribution by making the substitution (3) The corresponding cdf P(S) (fraction of amplitudes that are smaller than S) for the lognormal distribution is

WEIBULL DISTRIBUTION
The Weibull distribution is often used in engineering applications (Meeker and Escobar 1991;Weibull 1951). The pdf for the Weibull distribution is given by (Patel et al. 1976) where β and τ are fitting parameters. The mean μ and the aperiodicity (coefficient of variation) C V of the Weibull distribution are given by is the gamma function of x. The cdf for the Weibull distribution is given by If β = 1 the Weibull distribution becomes the exponential distribution with σ = μ and C V = 1. In the limit β → +∞ the Weibull distribution becomes exactly repetitive (a δ-function) with σ = C V = 0. In the range 0 < β <1 the Weibull distribution is often referred to as the stretched exponential distribution.

TESTS OF GOODNESS-OF-FIT
In order to determine whether a specific distribution is preferred, it is necessary to utilize tests of goodness-of-fit. Many such tests are available (Press et al. 1995). In this paper we quantify the goodness-of-fit of distributions using two tests.

KOLMOGOROV-SMIRNOV TEST
The first test we use is the Kolmogorov-Smirnov test (Press et al. 1995). To use this test the maximum absolute difference D KS between the cdf of the sample distribution (actual data) y i and the fitted distribution i y ) is determined: Then the significance level probability of the goodness-of-fit (the probability that the applied distribution is relevant) is given by and n is the number of data points. The preferred distribution has the smallest value of D KS and the largest value of Q KS .

ROOT MEAN SQUARED ERROR TEST
The second test of goodness-of-fit we use is the root mean squared error (RMSE). As it is suggested by the test name, it is the square root of the sum of squares of errors divided by the difference between the number of data points and the number of fitting parameters n is the number of data point, and k is the number of fitting parameters (k = 2 for our three distributions). This test is also known as the fit standard error and the standard error of the regression. The preferred distribution has the smallest RMSE value.

PARKFIELD SEQUENCE
Ideally, recurrent sequences of earthquakes should be used to establish a preferred statistical distribution. Unfortunately, the number of earthquakes available through historical records is too small for adequate statistical testing.
As an example we consider the sequence of seven characteristic earthquakes that occurred on the Parkfield (California) section of the San Andreas fault between 1857 and 2004 (Bakun et al. 2005). The slip rate is quite high ( ≈ 30 mm/year) and the earthquake magnitudes are relatively small (m ≈ 6.0), thus the recurrent times are short ( ≈ 25 years). Also, this fault is subject to a near constant tectonic drive due to the relative motion between the Pacific and North American plates. Slip on the Parkfield section of the San Andreas fault occurred in 1881, 1901, 1922, 1934, 1966, and 2004 with magnitudes from 6.0 to 6.05 by the instrumental estimate and from 5.9 to 6.1 by the estimate from MMI for an epicenter location on the 13 As the size of an event we can use here the seismic moment or energy of this event. However, the small number of registered earthquakes makes the application of statistical estimations impossible (Savage 1994). And this problem is relevant not only for the Parkfield sequence. Other earthquake sequences also have short statistics (e.g., Okada et al. 2003;Park and Mori 2007 an actual statistical distribution that should be used instead of the δ-function? This paper answers these questions.

ANDREAS FAULT
We now consider the recurrent statistics of slip events on the creeping section of the San Andreas fault, California. To do this we utilize records of two creepmeters on the San Andreas fault (Schulz 1989 show that the average long-term creep rate is in the range from 6 mm/year to 15 9 mm/year. The recorded data for both creepmeters can be downloaded from the U.S. Geological Survey web site (Langbein 2004).
For both creepmeters the data contain both daily and 10 minute telemetry records. Although the daily records are longer and give longer sequences, the 10 minute data also are used independently because they provide more accurate slip amplitude resolution.
Because the creepmeter piers are installed at shallow depths For earthquakes the duration of an event is of the order of seconds to minutes and can be neglected in comparison with the preceding period of slow stress accumulation. In contrast, for creep records, the duration of an event can be of the order of the interval between events. Therefore special techniques are required to distinguish one event from the next. To do this we 16 use the criterion of stationarity. If, after a well defined jump, the slip rate returns to a stationary creeping state before the next jump, these two jumps are treated as separate events. Otherwise, if one jump triggers another one in a transient, non-stationary process, these jumps are considered to be a single event.
Each creep record provides a unique opportunity to determine the complete sequence of events taking place at a given creepmeter location. In contrast to earthquakes, the rate of occurrence for creep events is much higher. Also each observed sequence of creep events provides the complete record of all events occurred at a given location. This gives an opportunity to associate events not only with the given fault but even with the given point of this fault. Therefore the reconstruction of event sequences from creep records gives both the longest possible sequences (up to 100 events) and the most accurate determination of recurrent events.
To filter the telemetry noise, threshold levels for slip amplitude are used: 0.077 mm for 10 minute xhr2 telemetry, 0.078 mm for 10 minute cwn1 telemetry, 0.3 mm for xhr2 daily telemetry, and 0.31 mm for cwn1 daily telemetry.
As the size of an event we will use the slip amplitude of this event.
The cumulative recurrent frequency-amplitude distributions of slip events for the 10 minute xhr2 and cwn1 records are given in Figs Finally, we integrate events from small to large amplitudes and divide the result by the total number of events.
The cumulative distribution of 51 recurrent slip amplitudes P(S) for the xhr2 10 minute event sequence is given as a function of the slip amplitude S in Fig. 2(a). The cumulative distribution of 45 recurrent slip amplitudes P(S) for the cwn1 10 minute event sequence is given as a function of the slip amplitude S in Fig. 2(c). The means and coefficients of 19 variation of these recurrent slip amplitudes are given in Table 1. Also included in Figs. 2(a) and 2(c) are the best-fits (maximum likelihood fits) of the Brownian passage-time, lognormal, and Weibull distributions. Both the parameters of these fits and the goodness-of-fit estimators are given in Table 1.
In Figs. 2(b) and 2(d) the recurrent statistics for the xhr2 and cwn1 10 minute event sequences are plotted in the form -ln(1-P(S)) versus S in log 10 -log 10 axes. In this form the Weibull distribution is a straight-line fit with slope β so that this is known as a Weibull plot. Also included in these Figures are the previous best-fits of the Brownian passage-time, lognormal, and Weibull distributions with RMSE in the log 10 -log 10 axes given in Table 1.
The cumulative distribution of 76 recurrent slip amplitudes P(S) for the xhr2 daily event sequence is given as a function of the slip amplitude S in Fig. 3(a). The cumulative distribution of 104 recurrent slip amplitudes P(S) for the cwn1 daily event sequence is given as a function of the slip amplitude S in Fig. 3(c). The means and coefficients of variation of these recurrent slip amplitudes are given in

SLIDER-BLOCK MODEL
In this section we consider the behavior of a slider-block model in order to study the statistics of event sizes. We utilize a variation of the linear sliderblock model which Carlson and Langer (1989) used to illustrate the selforganization of such models. We consider a linear chain of 500 slider blocks of mass m pulled over a surface at a constant velocity V L by a loader plate as illustrated in Fig. 4. Each block is connected to the loader plate by a spring with spring constant k L . Adjacent blocks are connected to each other by 21 springs with spring constant k C . Boundary conditions are assumed to be periodic: the last block is connected to the first one.
The blocks interact with the surface through friction. In this paper we prescribe a static-dynamic friction law. The static stability of each sliderblock is given by where F Si is the maximum static friction force on block i holding it motionless, and y i is the position of block i relative to the loader plate.
When the cumulative force of the springs connecting to block i exceeds the maximum static friction F Si , the block begins to slide. We include inertia, and the dynamic slip of block i is controlled by the equation where F Di is the dynamic (sliding) frictional force on block i. The loader plate velocity is assumed to be much smaller than the slip velocity, requiring It is convenient to introduce the non-dimensional variables and The ratio of static to dynamic friction φ is assumed to be the same for all In terms of these non-dimensional variables the static stability condition (13) becomes Before obtaining solutions, it is necessary to prescribe the parameters φ , α, and β i . The parameter α is a tuning parameter and is the stiffness of the system. We consider α = 1000 which corresponds to a very stiff model. We use this high value for the stiffness because the stiff slider-block model 24 obeys the principles of SOC. The ratio φ of static friction to dynamic friction is taken to be the same for all blocks 5 . 1 = φ , while the values of frictional parameters β i are assigned to blocks by uniform random distribution from the range 1 < β i < 3.5. This random variability in the system is a "noise" required to thermolize the system and generate event variability.
The loader plate springs of all blocks extend according to Eq. (19) until a block becomes unstable from Eq. (18). The dynamic slip of that block is calculated using the Runge-Kutta numerical method to obtain a solution of Eq. (20). A coupled 4 th -order iterational scheme is used, and all equations are solved simultaneously (the Runge-Kutta coefficients of neighboring blocks participate in the generation of the next order Runge-Kutta coefficient for the given block). The dynamic slip of one block may trigger the slip of other blocks and the slip of all blocks is followed until they all become stable. Then the procedure repeats. For the stiff system with α = 1000 the system-wide (500 block) events dominate. If we assume that each block represents an asperity and the whole model represents a fault or fault segment it is reasonable to conclude that the system-wide events correspond to characteristic events. First, we will consider the recurrent statistics at a given point on a fault. In the case of the slider-block model this corresponds to a given block of the model. We 25 consider statistics for the 'strongest', 'weakest', and 'medium' blocks. As a strongest block we choose the block with the highest coefficient of friction, i.e., a block with the highest β i . As a weakest block we choose the block with the lowest coefficient of friction. And as a medium block we choose the block with the friction coefficient which is close to the friction coefficient averaged over the model. As the size of an event we choose the total slip amplitude of the given block during this event. 1(c), and 1(d). We see here the same anomaly due to the presence of noncharacteristic (non-system-wide) events. However, now for the case of numerical simulations we have an opportunity to separate characteristic events directly as system-wide events (without the introduction of the amplitude threshold). Similar to the case of creep events, we remove noncharacteristic events from statistics, integrate cdf from small to large amplitudes, and normalize the statistics dividing it by the total number of events. 26 The cumulative distribution of recurrent slip amplitudes P(S) for the event sequences of the strongest, weakest, and medium blocks are given as functions of the slip amplitude S in Figs. 6(a), 6(c), and 6(e). The means and coefficients of variation of these sequences are given in Table 2. Also included in Figs. 6(a), 6(c), and 6(e) are the best-fits (maximum likelihood) of the Brownian passage-time, lognormal, and Weibull distributions. Both the parameters of these fits and the goodness-of-fit estimators are given in Table 2.
Figs. 6(b), 6(d), and 6(f) present the Weibull plots corresponding to Figs. 6(a), 6(c), and 6(e) respectively. Also included in these Figures are the corresponding best-fits of the Brownian passage-time, lognormal, and Weibull distributions with the RMSE in log 10 -log 10 axes given in Table 2.
The sequences investigated above have been the sequences at a given point of a fault. But for the slider-block model we can obtain sequences at a given fault. We consider again the system-wide events as the characteristic events. Then we can consider energy dissipated by the whole model (by all blocks) during an event as the size of this event. Indeed, the energy dissipated by all blocks during an event is already not associated with the slip amplitude at a given point of the model but is associated with the slip amplitude averaged over the model. The frequency-size statistics could be constructed both for energies and slip amplitudes -there are no scientific 27 objections to do it anyway. Therefore as the size of an event we can use the energy of this event as well as the slip amplitude.
The cumulative distribution of 715 recurrent energies P(E) for the sequence of system-wide events is given as a function of the dissipated energy E in Fig. 7(a). The mean and coefficient of variation of these recurrent energies are given in Table 3. Also included in Fig. 7(a) are the best-fits (maximum likelihood) of the Brownian passage-time, lognormal, and Weibull distributions. Both the parameters of these fits and the goodness-of-fit estimators are given in Table 3. Fig. 7(b) presents the Weibull plot corresponding to Fig. 7(a). Also included in this Figure are the corresponding best-fits of the Brownian passage-time, lognormal, and Weibull distributions with the RMSE in log 10log 10 axes given in Table 3.
The estimators of goodness-of-fits given above demonstrate convincingly that the fits of the Weibull distribution for all four sequences are much better than the lognormal or Brownian passage-time distributions. . Therefore, the creep events and the slider-block model exhibit closely related recurrent behavior of both frequency-size and time-interval statistics. Two independent systems whose applicability to earthquakes is often discussed in the literature have very closely related behavior. This gives us a hope that the stiff slider-block model, which obeys the principles of SOC, could represent the behavior of actual sequences of characteristic earthquakes at a given fault.

SAND-PILE MODEL
The stiff slider-block model obeys the principles of SOC. Applicability of this concept to earthquakes has been discussed in the literature. Therefore it would be interesting to take a look at the recurrent frequency-size behavior of a sand-pile model as a classical representative of SOC (Bak et al. 1988).
We utilize the simplest variation of 2-dimensional sand-pile model. A square grid of 100 by 100 sites is strewn by sand grains. When a site accumulates four or more grains it becomes unstable and redistributes four grains to its four neighbors. Instability of one site can trigger instability of others forming a complex avalanche in the system. The sand redistribution during an avalanche is assumed to be much faster than the rate of the stable sand accumulation during strewing (similarly, velocity of earthquake propagation is much faster than tectonic rate of stress accumulation). Therefore the slow sand strewing is neglected during avalanches. Boundaries of the model are 31 assumed to be free so the sand can leave the model at boundaries when it is redistributed by a boundary site.
To construct the recurrent statistics at a given point of the model we choose site in the middle of the lattice. All avalanches with participation of this site will be counted as events at this site of the model. To separate characteristic events we need to construct a criterion similar to the systemwide criterion for the slider-block. For the sand-pile model we choose a percolation criterion as a criterion for an event to be characteristic. I.e., an event will be considered as characteristic if it percolates the lattice and connects all four boundaries (up-right-down-left percolation). As the size of an event we choose the number of different sites participating in this avalanche. Each site can repeatedly loose its stability during this event but will be counted only once in the size of the event.
The cumulative distribution of 2013 recurrent events P(S) is given as a function of the event size S in Fig. 8(a). The mean and coefficient of variation of this statistics are given in Table 4. Also included in Fig. 8(a) are the best-fits (maximum likelihood) of the Brownian passage-time, lognormal, and Weibull distributions. Both the parameters of these fits and the goodness-of-fit estimators are given in Table 4.  Table 4.
Again, the Weibull distribution is the best-fit distribution. However, now its exponent β has much higher value β = 8.09 which corresponds to the aperiodicity C V = 0.15. This is probably the effect of the percolation criterion as a criterion for an event to be characteristic. Although for this case the concept of self-organized criticality preserves the same functional (Weibull) dependence of the recurrent frequency-size statistics, the sand-pile model has another symmetry and belongs to another universality class with different anomalous dimensions. This question requires further investigation. . The same tendency was also found by Abaimov et al. (2007a;2007b) for the recurrent time-interval behavior of these two systems. The fact that two independent systems whose applicability to earthquakes is often discussed in the literature have closely related recurrent behavior supports the point of view that these two systems actually represent the recurrent earthquake behavior at a given fault or fault segment.   Weibull plot is given in (b) as diamonds. In both cases the data are compared with the best-fit Brownian passage-time distribution (dash-dot lines), the best-fit lognormal distribution (short-dash lines), and the best-fit Weibull distribution (long-dash lines).