Apr 24, 2014 - The upper tail of the Îº-Weibull distribution declines as a power law in contrast with. Weibull scaling. ...... and metadata required i...

0 downloads 0 Views 887KB Size

arXiv:1308.1881v3 [physics.geo-ph] 24 Apr 2014

Technical University of Crete, Chania 73100, Greece Giorgio Kaniadakis‡ Department of Applied Science and Technology, Politecnico di Torino, Corso Duca degli Abruzzi 24, 10129 Torino, Italy (Dated: April 25, 2014)

1

Abstract The Weibull distribution is a commonly used model for the strength of brittle materials and earthquake return intervals. Deviations from Weibull scaling, however, have been observed in earthquake return intervals and the fracture strength of quasi-brittle materials. We investigate weakest-link scaling in finite-size systems and deviations of empirical return interval distributions from the Weibull distribution function. Our analysis employs the ansatz that the survival probability function of a system with complex interactions among its units can be expressed as the product of the survival probability functions for an ensemble of representative volume elements (RVEs). We show that if the system comprises a finite number of RVEs, it obeys the κ-Weibull distribution. The upper tail of the κ-Weibull distribution declines as a power law in contrast with Weibull scaling. The hazard rate function of the κ-Weibull distribution decreases linearly after a waiting time τc ∝ n1/m , where m is the Weibull modulus and n is the system size in terms of representative volume elements. We conduct statistical analysis of experimental data and simulations which shows that the κ-Weibull provides competitive fits to the return interval distributions of seismic data and of avalanches in a fiber bundle model. In conclusion, using theoretical and statistical analysis of real and simulated data, we demonstrate that the κ-Weibull distribution is a useful model for extreme-event return intervals in finite-size systems. PACS numbers: 02.50.-r, 89.75.Da, 62.20.mj Keywords: waiting times, heavy tails, extreme events, Weibull

∗

[email protected]

†

[email protected]

‡

[email protected]

2

I.

INTRODUCTION

Extreme events correspond to excursions of a random process X(t), where t is the time index, to values above or below a specified threshold zq . In natural processes, extreme events include unusual weather patterns, ocean waves, droughts, flash flooding, and earthquakes. Such phenomena have important social, economic and ecological consequences. The FisherTippet-Gnedenko (FTG) theorem states that if {Xi }ni=1 are independent and identically distributed (i.i.d.) variables, then a properly scaled affine transformation of the minimum χn := min(X1 , . . . , Xn ) follows asymptotically (for n → ∞) one of the extreme value distributions, which include the Gumbel (infinite support), reverse Weibull [1] (positive support) and Fr´echet distributions (negative support) [2]. Whereas the FTG theorem is a valuable starting point, many processes of interest involve complex systems with correlated random variables. The impact of correlations on the statistical behavior of complex physical systems thus needs to be understood. Early research on extreme events statistics focused on purely statistical approaches [3, 4]. Current efforts are based on nonlinear stochastic models and aim to understand the patterns exhibited by extreme events and to control them [5–10]. To improve risk assessment methodologies, the statistics of the return intervals, i.e., the time that elapses between consecutive crossings of a given threshold by X(t), is an important property. If the threshold crossing implies failure (e.g., fracture), then the return intervals are intimately linked to the strength distribution of the system [11]. Herein we focus on the return intervals of earthquakes, i.e., earthquake return intervals (ERI) [12] and the return intervals of avalanches in fiber bundle models under compressive loading. From a broader perspective, our scaling analysis can be also applied to other systems or properties governed by weakest-link scaling laws, such as the strength of quasibrittle heterogeneous materials. This document is structured as follows. In the remainder of this section we review the literature on earthquake return intervals. Section II presents the basic principles of weakestlink scaling and its connection to the Weibull distribution. In Section III we present an extension of weakest-link scaling for finite-size systems and motivate the use of the κ-Weibull distribution. Section IV links the κ-Weibull distribution to earthquake return intervals using theoretical arguments. In Section V we apply these ideas to seismic data. Section VI focuses on the return intervals between avalanches in a fiber bundle model with global load sharing and demonstrates the performance of the κ-Weibull distribution on this synthetic data. 3

Finally, Section VII summarizes our conclusions and briefly discusses the significance of the results.

A.

Earthquake Return Intervals

Earthquake patterns can be investigated over different spatial supports which range from a single fault to a system of faults [13]. Both isolated faults and fault systems represent complex problems that combine nonlinear and stochastic elements. Various probability functions have been proposed to model the earthquake return interval distribution (for a recent review see [14]). Several authors have proposed that earthquakes are manifestations of a self-organized system near a critical point [15–17] or of a system near a spinodal critical point [18, 19]. Both cases imply the emergence of power laws. Bak et al. [15] introduced a global scaling law that relates earthquake return intervals with the magnitude and the distance between the earthquake locations. These authors analyzed seismic catalogue data over a period of 16 years from an extended area in California that includes several faults (ca. 3.35 × 105 events). They observed power-law dependence over eight orders of magnitude, indicating correlations over a wide range of return intervals, distances and magnitudes. Corral and coworkers [16, 20–22] introduced a local modification of the scaling law so that the return intervals probability density function (pdf) follows the universal expression fτ (τ ) ' λf˜(λ τ ), where f˜(τ ) is a scaling function and the typical return interval τ¯ is specific to the region of interest. Saichev and Sornette [17, 23] generalized the scaling function by incorporating parameters with local dependence. Their analysis was based on the mean-field approximation of the return intervals pdf in the epidemic-type aftershock sequence (ETAS) model [24]. ETAS incorporates the main empirical laws of seismicity, such as the Gutenberg-Richter dependence of earthquake frequency on magnitude, the Omori-Utsu law for the rate of the aftershocks, and a similarity assumption that does not distinguish between foreshocks, main events and aftershocks (any event can be considered as a trigger for subsequent events). Several studies of earthquake catalogues and simulations show that the Weibull distribution is a good match for the empirical return intervals distribution [25–34]. In addition to statistical analysis, arguments supporting the Weibull distribution are based on Extreme Value Theory [35], numerical simulations of slider-block models [32], and growth-decay mod4

els governed by the geometric Langevin equation [6]. The Weibull distribution is also used to model the fracture strength of brittle and quasibrittle engineered materials [36–38] and geologic media [39]. With respect to extreme value theory, if we ignore correlations the FTG theorem favors the Weibull because the return intervals are non-negative, whereas the Fr´echet distribution for minima has negative support and the Gumbel distribution has unbounded support. A physical connection between the distribution of shear strength of the Earth’s crust and the ERI distribution was proposed in [11]. According to a simplified stick-slip model, if the shear strength follows the Weibull distribution, under certain conditions the ERI also follows the Weibull distribution with parameters which are determined from the respective strength parameters and the exponent of the loading function. The conditions include: (i) the stress increase during the stick phase follows a power-law function of time (ii) the duration of the slip phase can be ignored (iii) the residual stress is uniform across different stick-slip cycles, and (iv) the parameters of the Earth’s crust shear strength distribution are uniform over the study area. In particular, if the shear strength follows the Weibull distribution with modulus ms and the stress increases with time as a power law with exponent β between consecutive events, then the ERIs also follow the Weibull distribution with modulus m = ms β. On a similar track, a recent publication reports strong connections between the statistics of laboratory mechanical fracture experiments and earthquakes [40, 41].

II.

WEAKEST-LINK SCALING

The weakest-link scaling theory underlies the Weibull distribution. Weakest-link scaling was founded by the works of Gumbel [3] and Weibull [4] on the statistics of extreme values; it is used to model the strength statistics of various disordered materials [36, 42–44]. Weakestlink scaling treats a disordered system as a chain of critical clusters, also known as links or representative volume elements (RVEs). The strength of the system is determined by the strength of the weakest link, hence the term weakest-link scaling [45]. The concept of links is straightforward in simple systems, such as one-dimensional chains. In higher dimensions the RVEs correspond to critical subsystems, possibly with their own internal structure, failure of which destabilizes the entire system [46]. We consider systems that follow weakest-link scaling and comprise n links. We use the symbol x to denote the values of a random variable 5

X which can represent mechanical strength or time intervals between two events. (i)

We denote by F1 (x) = Prob(X ≤ x) the cumulative distribution function (cdf ) that X takes values that do not exceed x. For example, if X denotes mechanical strength (return (i)

intervals), then F1 (x) is the probability that the i-th link has failed when the loading has reached the value x (when time τ = x has passed). Respectively, we denote by Fn (x) the probability that the entire system fails at x. The function Rn (x) := 1 − Fn (x) represents the system’s survival probability. The principle of weakest-link scaling is equivalent to the statement that the system’s survival probability is equal to the product of the link survival probabilities; this is expressed mathematically as Rn (x) =

n Y

(i)

R1 (x).

(1)

i=1 (i)

If all the RVEs share the same functional form for R1 (x), (1) leads to Rn (x) =

h

(i) R1 (x)

in

.

(2)

(i)

Assuming that R1 (x) is independent of n, Eq. (2) implies the following scaling expression for n > n0 0

Rn (x) = [Rn0 (x)]n/n .

(3)

(i)

If the Weibull ansatz R1 (x) = e−φ(x) is satisfied [4], then Rn (x) = e−n φ(x) . Furthermore, if φ(x) = (x/x0 )m , then F (x) := F (i) (x) for i = 1, . . . , n, and x

m

F (x) = 1 − e−( xs ) ,

(4)

where xs is the scale parameter and m > 0 is the Weibull modulus or shape parameter. The size dependence of xs is determined by xs = x0 /n1/m [47]. Let us define the double logarithm of the inverse of the survival function Φn (x) = ln ln Rn−1 (x). In light of (3), the following size-dependent scaling is obtained Φn (x) = Φn0 (x) + ln(n/n0 ).

(5)

Based on the weakest-link scaling relation (5) and the pioneering works [48, 49], it can be shown using asymptotic analysis that the system’s cdf tends asymptotically (as n → ∞) to the Weibull cdf [49, 50]. Curtin then showed that the large-scale cdf parameters depend both on the system and the RVE size [42]. 6

The Weibull pdf is given by f (x) = dF (x)/dx and leads to the expression m−1 x m x e−( xs ) . f (x) = m τs

(6)

For m < 1 the Weibull is also known as the stretched exponential distribution [5] and finds applications in generalized relaxation models [51, 52], whereas for m = 2 it is equivalent to the Rayleigh distribution. For m < 1 the pdf has an integrable divergence at x = 0 and decays exponentially as x → ∞. For m = 1 the exponential pdf is obtained, whereas for m > 1 the pdf develops a single peak with diminishing width as m ↑. Finally, for the Weibull distribution the function Φn (x) is linearly related to the logarithm 1/m 0 . of x, i.e., Φn (x) = m ln(x/xs ), and (5) implies the size dependence xxss = nn0

III.

WEAKEST-LINK SCALING AND FINITE-SIZE SYSTEMS

The Weibull model assumes the existence of independent RVEs and n 1. Nevertheless, there are systems for which the asymptotic assumption n 1 is not a priori justified. For example, fault systems span a wide range of scales (100 - 106 m). The size or even the existence of an RVE are not established for fault systems. In quasibrittle materials, the RVE is assumed to exist but its size is not negligible compared to the system size, leading to deviations from the Weibull scaling in the upper tail of the strength pdf [37, 38, 46, 53]. Using a piecewise Weibull-Gaussian model for the RVE strength pdf, Bazant et al. [37, 38, 46, 53] proposed that the system pdf exhibits a transition from Weibull scaling in the lower (left) tail to Gaussian dependence in the upper tail at a probability threshold that moves upward as the size increases. We consider a system that follows weakest-link scaling and consists of RVEs with uniform properties. We associate the parameter κ with the number of effective RVEs through n = 1/κ. Hence, κ (and also n) are parameters to be estimated from the data. Note that n does not need to be integer, whereas for systems smaller than one RVE n < 1 (κ > 1) is possible.

A.

κ-Weibull Distribution

The exponential tail of the Weibull pdf defined in (6) follows from the fact that the survival probability R(x) = exp(−[x/xs ]m ) is defined in terms of the exponential function 7

exp(−z). On the other hand, in the last decades particular attention has been devoted to pdfs that exhibit power-law tails, namely Ax−α . Such dependence has been observed in many branches of natural sciences including seismology, meteorology, and geophysics [54, 55]. The simplest way to treat systems with these features is to replace the exponential function in the definition of R(x) by another proper function which generalizes the exponential function and presents power-law tails. A one-parameter generalization of the exponential function has been proposed in [56, 57] and is given by 1/κ √ 1 + z 2 κ2 + zκ , expκ (z) =

(7)

with 0 ≤ κ < 1. The above generalization of the ordinary exponential emerges naturally within the framework of special relativity, where the parameter κ is proportional to the reciprocal of light speed [58, 59]. In that context, expκ (z) is the relativistic generalization of the classical exponential exp(z). The inverse function of the κ-exponential is the κ-logarithm, defined by lnκ (z) =

z κ − z −κ . 2κ

(8)

By direct inspection of the first few terms of the Taylor expansion of expκ (z), reported in [60] expκ (z) = 1 + z +

z2 z3 z4 + (1 − κ2 ) + (1 − 4κ2 ) + . . . , 2 3! 4!

(9)

it follows that when z → 0 or κ → 0 the function expκ (z) approaches the ordinary exponential i.e. expκ (z) ∼ exp(z),

(10a)

expκ (z) ∼ exp(z).

(10b)

z→0

κ→0

The most important feature of expκ (z) regards its power-law asymptotic behavior [57, 60] i.e. ±1/κ expκ (z) ∼ 2κz .

(11)

z→±∞

We remark that the function expκ (−z) for z → 0 coincides with the ordinary exponential i.e. expκ (−z) ∼ exp(−z), whereas for z → +∞ it exhibits heavy tails i.e. expκ (−z) ∼ (2κz) −1/κ . Therefore the function expκ (−z) is particularly suitable to define the survival probability [61, 62]. Following the change of variables κ = 1/n and z = (x/xs )m we obtain Rκ (x) = expκ (− [x/xs ]m ) , 8

(12)

0

10

−1

10

κ=0 κ=0.1 κ=0.5 κ=0.9 κ=2 κ=3

0.1

0.08

−2

f(x)

f(x)

0.12

κ=0 κ=0.1 κ=0.5 κ=0.9 κ=2 κ=3

10

0.06

0.04 −3

10

0.02 −4

10

0

20

40

60

80

100

0 0

5

10

15

x

(a) m = 0.7

20 x

25

30

35

40

(b) m = 3

FIG. 1: κ-Weibull pdfs for xs = 10, different values of κ, and (a) m = 0.7 (b) m = 3. The resulting κ-Weibull distribution exhibits a power-law tail inherited by the κ-exponential: Fκ (x) = 1 − expκ (− [x/xs ]m ) , m−1 m x exp (− [x/xs ]m ) p κ . fκ (x) = xs xs 1 + κ2 (x/xs )2m

(13) (14)

Plots of the κ-Weibull pdf for xs = 10, different κ, and two values of m (m < 1 and m > 1) are shown in Fig. 1. The plots also include the Weibull pdf (κ = 0) for comparison. For both m higher κ lead to a heavier right tail. For m = 3 the mode of the pdf moves to the left of xs as κ increases. To the right of the mode, lower κ correspond, at first, to higher pdf values. This is reversed at a crossover point beyond which the higher-κ pdfs exhibit slower power-law decay for x → ∞, i.e., fκ (x) ∝ x−α , where α = 1 + m/κ. The crossover point occurs at ≈ 1.5xs for m = 3, whereas for m = 0.7 at ≈ 5 xs . For m = 0.7 the mode is at zero independently of κ, since the distribution is zero-modal for m ≤ 1. It is important to note that the κ-Weibull admits explicit expressions for all the important univariate probability functions. The κ-Weibull hazard rate function is defined by means of hκ (x) = fκ (x)/Rκ (x) = −d ln Rκ (x)/dx, leading to hκ (x) =

m (x/xs )m−1 p . xs 1 + κ2 (x/xs )2m

The κ-Weibull quantile function for a given survival probability r is defined by 1/m 1 1 Tκ (r) = lnκ . xs r 9

(15)

(16)

In addition, if we define Φ0 κ (x) = ln lnκ (1/Rκ (x)) it follows that Φ0 κ (x) = m ln

x xs

. Hence,

Φ0 κ (x) is independent of κ and regains the logarithmic scaling of the double logarithm of the inverse survival function.

B.

RVE Survival Function

We define the RVE cdf at level x ∈ [0, ∞) through the equation m s 2m 1 x x 1 F1 (x) = 1 + − 1+ 2 . n xs n xs

(17)

F1 (x) is a well defined cdf, because F1 (x = 0) = 0, whereas for x > 0 F1 (x) is an increasing function of x, and limx→∞ F1 (x) = 1. This particular form of F1 (x) is motivated by arguments similar to those used in the Weibull case. In Section II, the Weibull survival function Rn (x) (i)

was derived from (1) assuming that the link survival function is R1 (x) = e−φ(x) where φ(x) = (x/xs )m . Another approach that does not require the exponential dependence of the RVE survival function is based on the following approximation Rn (x) = [1 − F1 (x)]n ⇒ ln Rn (x) = n ln [1 − F1 (x)] ≈ −n F1 (x). The above assumes that F1 (x) 1 for the link cdfs if n is large. Then, assuming that F1 (x) ∝ (x/xs )m the Weibull form is obtained. The dependence of F1 (x) for large x which becomes relevant for finite n, however, is not specified. In contrast, (17) generalizes the algebraic dependence so that F1 (x) ∼ (x/xs )m for x → 0, whereas F1 (x) is also well defined for x → ∞. From (17) it follows that the respective survival function is s 2m m 1 x 1 x R1 (x) = 1 + 2 − . n xs n xs

(18)

Application of the weakest-link scaling relation (2) to (18) leads to the following system survival function

s Rn (x) =

1+

1 n2

x xs

2m −

1 n

x xs

m

n .

(19)

The definition (17) implies that F1 (x) and R1 (x) depend on the number of RVEs, which destroys the weakest-link scaling relation (3). Based on (18) and using z = (x/xs )m it follows 10

that ∂R1 (z) z = 2 ∂n n

z 1− √ > 0, ∀z > 0. z 2 + n2

Hence, the survival probability of single RVEs at a given threshold z increases with n. We propose an ansatz which is consistent with the dependence of R1 (x) as given by (18). Assume that the system comprises a number of units (e.g., faults) with inter-dependent RVE survival probabilities, as expected in the presence of correlations. Following a renormalization group (RG) procedure, the interacting units are replaced by non-interacting “effective RVEs”. The RG procedure recovers the product form (2) for the survival probability of independent RVEs, while renormalizing the scale parameter x0 by the number of effective RVEs. We can think of κ = 1/n as a measure of the range of interactions versus the size of the system; κ = 0 yields the classical Weibull pdf for infinite systems, whereas κ ↑ implies that the range of correlations increases thus reducing the number of independent units; the case κ = 1 means that the system can not be reduced to smaller independent units.

IV.

WEAKEST-LINK SCALING AND RETURN INTERVALS

Below we focus explicitly on earthquake return intervals; thus, we replace x with τ . In earthquake analysis the spatial support includes either a single fault or a system of several faults. The notion of an RVE with respect to earthquakes is neither theoretically developed nor experimentally validated. Hence, herein we assume that the study domain involves n independent, identically distributed RVEs, where n is not necessarily an integer [64]. An earthquake catalog is a table of the marked point process [65] C = {si , ti , Mi }i=1,...,N , where si is the location, ti the time, and Mi the magnitude of the seismic event. Given a threshold magnitude Mc , an ERI sequence comprises the intervals {τj = tj+k − tj : (Mj , Mj+k > Mc ) ∧ (Mj+1 , . . . , Mj+k−1 ≤ Mc )}, where j = 1, . . . , Nc − 1, Nc is the number of events with magnitude exceeding Mc (Fig. 2), and ∧ is the logical conjunction symbol. (i)

The random variable TMc (i = 1, . . . , n) denotes the quiescent interval for the i-th RVE (i)

during which no events of magnitude M > Mc occur. The cdf F1 (τ ; Mc ) = Prob(TMc ≤ τ ) represents the probability of RVE “failure”, i.e., that an event with M > Mc occurs on the RVE within time interval τ from the previous event. In the following, we suppress the dependence on Mc for brevity. 11

M Mj

Mj+k τj

Mc

tj

tj+k

t

FIG. 2: Schematic illustrating the definition of a return interval τj as the time between two consecutive events with magnitudes Mj and Mj+k that exceed the threshold Mc , i.e., the events occurring at times tj and tj+k . A.

Survival Probability Function

The survival probability R1 (τ ) = 1−F1 (τ ) is the probability that no event with magnitude (i)

M > Mc occurs on the RVE during the interval TMc ≤ τ . For τ → 0, it follows from (18) that R1 (τ ) ∼ 1 − τ m /(n τsm ). For τ 1 and finite n, it follows that R1 (τ )∼n τsm /2τ m , and thus Rn (τ ) shows the power-law dependence Rn (τ )∼(n/2)n (τs /τ )mn , characteristic of the Pareto distribution. In addition, limn→∞ Rn (τ ) = exp [−(τ /τs )m ], thus recovering the Weibull survival probability at the limit of an infinite system. The above equation shows that the interval scale for large n saturates at τs , in contrast with the classical τs ∝ n−1/m Weibull scaling. Based on the Gutenberg-Richter law of seismicity which predicts exponential decay of earthquake events as Mc ↑, it follows that τs ↑ as Mc ↑. In contrast, m is expected to vary more slowly with Mc [11].

B.

Median of Return Intervals

The median of the single RVE distribution is defined by R1 (τmed;1 ) = 0.5, and based 1/m on (18) it is given by τmed;1 = τs 3n . The median of the κ-Weibull distribution [63] 4 for a system of n = 1/κ RVEs is given by τmed;n = (lnκ 2)1/m τs , whereas the median of the Weibull distribution is lim τmed;n = (ln 2)1/m τs . Based on the above, the ratio of the median n→∞

return interval for a finite system over the median return interval of an infinite system both of which have the same τs , is given by τmed;n /τmed;∞ = (lnκ 2/ ln 2)1/m . The ratio is plotted in Fig. 3. For n fixed the ratio is reduced with increasing m, whereas for m > 1 the median 12

τmed;n/τmed;inf

2 1.8 1.6 1.4 1.2

2 4

8

10

m

6

4

2

n

FIG. 3: Median ratio τmed;n /τmed;∞ for the κ-Weibull distribution versus the Weibull modulus m and the system size n. return interval varies only slightly with m. Keeping m fixed, the median return interval ratio declines with n toward 1. This means that smaller systems have higher median return interval than the infinite system —assuming that the characteristic interval does not change with size. This result is related to the heavier (i.e., power-law) upper tail of the finite-size system.

C.

Hazard Rate Function

A significant question for seismic risk assessment is whether the probability of an earthquake of given magnitude grows or declines as the waiting time increases [66, 67]. An answer to this question involves the hazard rate function of the return intervals. The latter is the 13

conditional probability that an earthquake will occur at time τ ∗ within the infinitesimal time window τ < τ ∗ ≤ τ + dτ , given that there are no earthquakes in the interval [0, τ ]. Hence [66], hτ (τ ) =

fτ (τ ) Prob[τ < τ ∗ ≤ τ + dτ |τ ∗ > τ ] = . dτ R(τ )

If earthquakes were random (memoryless) processes, distributed in time according to the Poisson law, the ERI would follow the exponential distribution leading to a constant hτ (τ ). If the ERI follows the Weibull distribution with cdf (4), the hazard rate is given by m−1 τ m hτ (τ ) = . τs τs

(20)

According to (20), the hazard rate for m > 1 increases as τ → ∞. This is believed to apply to characteristic earthquakes that occur on faults located near plate boundaries. In contrast, the Weibull distribution with m < 1 as well as the lognormal and the power-law distributions exhibit the opposite trend [66]. Since Bak proposed a connection between earthquakes and self-organized criticality [15], universal or locally modified power-law expressions and the gamma probability density function —which is a power law with an exponential cutoff for large times— have been proposed as models of the ERI pdf [21, 23, 40, 68]. The behavior of the gamma distribution depends on the value of the power-law exponent in the same way as the Weibull model. An analysis of two earthquake catalogues based on the gamma distribution concludes that the hazard rate decreases with time (corresponding to an exponent between 0 and 1) [67]. The hazard rate of the κ-Weibull is given by (15). For finite n and for τ τs n1/m , hτ (τ ) ∼ 1/τ . If we take the limit n → ∞ before τ → ∞, the Weibull hazard rate (20) is obtained. For a fixed RVE size, even if n 1, the Weibull scaling holds for τ < τs n1/m whereas for τ τs n1/m the τ −1 scaling dominates. This behavior of h(τ ) is demonstrated in Fig. 4: For m < 1 the dependence is not severely affected by size effects; for m = 1 there is a constant plateau followed by an 1/τ decay, whereas for m > 1 the initial increase of h(τ ) turns into an 1/τ decay after a turning point which occurs for τc ≈ τs n1/m .

V.

ANALYSIS OF EARTHQUAKE RETURN INTERVAL DATA

The estimation of the ERI distribution from data is complicated by the fact that the κWeibull distribution and the Weibull distribution are close over the range 0 ≤ τ ≤ τw , where 14

4

10

3

10

W: m = 5 K: m = 5 W: m = 2 K: m = 2 Exponential K: m = 1 W: m = 0.5 K: m = 0.5

2

10

1

h(τ)

10

0

10

−1

10

−2

10

−3

10

0

10

1

2

10

10 τ/τ 0

3

10

4

10

FIG. 4: Log-log plot of hazard rates h(τ ) versus τ for different models. (K): κ-Weibull hazard rate based on (15) for n = 100 with m = 0.5 (magenta solid line and diamonds), m = 1 (black solid line), m = 2 (purple solid line with squares), and m = 5 (cyan line with circles). (W): Weibull hazard rate obtained from (20) with m = 0.5 (magenta dashed line with diamonds), m = 2 (purple dashed line with squares), m = 5 (cyan dashed line with circles). (Exponential): Exponential hazard rate obtained from (20) with m = 1 (black dash-dot line).

τw is a parameter that depends on n and m. Differences in the tail of ERI distributions are best visualized in terms of Φn (τ ), as shown in the Weibull plots of Figs. 5-6. On these diagrams, the deviation of the κ-Weibull distribution from the straight line diminishes with increasing n. The gamma distribution is also included for comparison purposes. The gamma probability model with pdf f (τ ) = τ α−1 exp(−τ /b) bα Γ(α) is often used in studies of earthquake return intervals, e.g. [16],[67]. For m < 1 the Weibull plot of the gamma 15

4 3

KW: N=1 KW: N=10 KW: N=100 Gamma Weibull

2

Φ(τ)

1 0 −1 −2 −3 −4 −5 −8

−6

−4

−2

0

2

4

6

ln(τ/τs) FIG. 5: Φn (τ ) versus ln(τ /τs ) for the Weibull distribution with m = 0.7 (green solid line), the κ-Weibull with m = 0.7 and n = 1 (blue dashed line), n = 10 (red dashed and dotted line), n = 100 (black dotted line), and the gamma distribution with α = 0.7 (magenta dashed line with small circles).

probability distribution is a convex function, whereas for m > 1 it becomes concave. In contrast, the κ-Weibull distribution is concave for all m.

A.

Microseismic sequence from Crete

We consider the return intervals for an earthquake sequence from the island of Crete (Greece) which involves over 1 821 micro-earthquake events with magnitudes up to 4.5 (ML ) (Richter local magnitude scale) [11]. The sequence was recorded between July 2003 and June 2004 [11]. The return intervals between successive earthquake events range from 16

6 4 2

Φ(τ)

0 −2 −4

KW: N=1 KW: N=10 KW: N=100 Gamma Weibull

−6 −8 −10 −3

−2

−1

0 ln(τ/τs)

1

2

3

FIG. 6: Φn (τ ) versus ln(τ /τs ) for the Weibull distribution with m = 3 (green solid line), the κ-Weibull with m = 3 and n = 1 (blue dashed line), n = 10 (red dashed and dotted line), n = 100 (black dotted line), and the gamma distribution with α = 3 (magenta dashed line with small circles).

1 (sec) to 19.5 (days). The spatial domain covered is approximately between 24.5◦ – 27◦ (East longitude) and 34◦ – 35.5◦ (North latitude). The magnitude of completeness for this data set is around 2.2 – 2.3 (ML ), which means that all events exceeding this magnitude are registered by the measurement network. We use the method of maximum likelihood to estimate the parameters of test probability distributions for the return intervals. The optimal κ-Weibull distribution for events above Mc = 2.3 (ML ) is compared with the optimal Weibull distribution in Fig. 7 [69]. Note that the empirical distribution of the return intervals has m < 1 and a concave tail, in contrast with the gamma density model (cf. Fig. 5). The κ-Weibull distribution approximates better 17

1 2 0.8

Φ(τ)

0 Empirical kappa−Weibull Weibull

−2

0.4

−4

0.2

−6

−8 4

0.6

6

8

10 ln(τ)

12

14

0 16

FIG. 7: Φ(τ ) versus ln(τ ) for earthquake return intervals of the Cretan earthquake sequence (CES). The return intervals refer to 628 events with magnitudes exceeding 2.3 (ML ). The magnitude of completeness is ≈ 2.2 (ML ). The maximum likelihood estimates of the κ-Weibull parameters are τˆ0 ≈ 3.19 × 104 (sec), m ˆ ≈ 0.78, κ ˆ ≈ 0.33. The estimate of Φ(τ ) using the empirical (data-based) cdf is shown with the solid blue line. The optimal κ-Weibull fit is shown with the red dashed line, whereas the optimal Weibull fit is shown with the green dashed and dotted line. The left vertical axis measures Φ(τ ), whereas the right vertical axis marks the corresponding cdf values.

the upper tail of the return intervals than the Weibull distribution. Both the Weibull distribution and the κ-Weibull distribution have a lighter lower tail than the data. These trends persist as Mc ↑, but the differences between the distributions progressively decrease. Fig. 8 compares the quantiles of the data distribution with those of the optimal κ-Weibull model. We investigate different hypotheses for the ERI distribution using the KolmogorovSmirnov test following the methodology described in [70].

The Kolmogorov-Smirnov

distance between the empirical (data) distribution, Femp (τ ) and the estimated (model) distribution, Fˆ (τ ), is given by D = supτ ∈R |Femp (τ ) − Fˆ (τ )|, where supA f (τ ) denotes the 18

5

kappa−Weibull (sec)

10

4

10

3

10

2

10

2

10

3

10

4

10 Data ERT (sec)

5

10

FIG. 8: Sample quantiles versus the corresponding quantiles of the optimal κ-Weibull distribution for the Cretan earthquake sequence return intervals (CES ERI) —Φ(τ ) of the data shown in Fig. 7. The κ-Weibull quantiles are obtained from (16) using the maximum likelihood estimates of the κ-Weibull parameters (see caption of Fig. 7). supremum of f (τ ) for τ ∈ A. The parameters of Fˆ (τ ) are also estimated using the method of maximum likelihood as described above. The null hypothesis is that Fˆ (τ ) represents the probability distribution of the data. We apply the test to the Poisson, normal (Gaussian), lognormal, Weibull, κ-Weibull, gamma, and generalized gamma distributions. The generm alized gamma distribution [71], with pdf given by f (x) = (d/xs )m xd−1 e−(x/xs ) Γ(d/m), incorporates both the gamma distribution (for m = 1) and the Weibull distribution (for m = d). The Kolmogorov-Smirnov test for a probability model with estimated parameters should be applied using Monte Carlo simulation to generate synthetic data from the estimated probability model. We generate random numbers from the Poisson, normal, lognormal, Weibull, and gamma distributions with the respective MATLAB random number generators. For the κ-Weibull (κ > 0) and for the generalized gamma distribution we implemented the inverse transform sampling method (see Appendix A). For each realization of return inter19

vals, we estimate the parameters of the optimal distribution model, Fˆ (j) (τ ), j = 1, . . . , Nsim . The Kolmogorov-Smirnov distance between the empirical distribution of the specific real(j)

ization, Femp (τ ), and the estimated model distributions for the j th realization are given (j) by D(j) = supτ ∈R |Femp (τ ) − Fˆ (j) (τ )|, j = 1, . . . , Nsim . The p-value of the KolmogorovPNsim 1 (j) Smirnov distance is defined as p = Nsim i=j 1D(j) >D (D ) —where 1A (τ ) = 1, if τ ∈ A and 1A (τ ) = 0, if τ ∈ / A, is the indicator function of the set A. The p-value is the probability the Kolmogorov-Smirnov distance will exceed D purely by chance if the null hypothesis is true. If p > pcrit , the null hypothesis is accepted, otherwise it is rejected. If we focus on the return intervals between earthquakes with magnitudes exceeding 2.3, the Kolmogorov-Smirnov test (based on 1 000 simulations) rejects the normal, gamma, and lognormal and Poisson distributions at the 5% significance level. In contrast, the Weibull and κ-Weibull distributions are accepted with p ≈ 0.09 and p ≈ 0.75, respectively, whereas the generalized gamma at p ≈ 0.10. For 3.5 ≥ Mc ≥ 2.3 (ML ) the calculated p-values are shown in Table I. The following remarks summarize the tabulated results: (i) The κ-Weibull and the Weibull distributions are accepted for all Mc except for 2.7 and 3.5. (ii) The gamma distribution is accepted for Mc ≥ 2.5 whereas the generalized gamma is accepted for all Mc except for Mc = 2.7. (iii) The lognormal is marginally accepted for Mc = 3.5. (iv) The Poisson distribution is accepted for Mc = 3.7, 3.9 (v) The normal and Poisson models are rejected for all Mc ≤ 3.7. (vi) For Mc ≤ 3.3 the Weibull, the gamma, generalized gamma and the κ-Weibull models have the highest p-values but their relative ranking changes with Mc . Note that for Mc > 3.3 the sample size is quite small (45 ≥ Nc ≥ 18) thus prohibiting the observation of long tails. Since for most Mc more than one model hypotheses pass the Kolmogorov-Smirnov test, it is desirable to somehow compare the different probability models. For this purpose we use the Akaike Information Criterion (AIC) [72]. The AIC is defined by AIC = 2 NLL + 2k, where NLL is the negative log-likelihood of the data for the given model, and k is the number of model parameters (k = 1 for the Poisson, k = 2 for the normal, lognormal, Weibull, and gamma, whereas k = 3 for the κ-Weibull and the generalized gamma). The term 2k in AIC penalizes models with more parameters. In general, a model with lower AIC is preferable to one with higher AIC. We present AIC results for the Crete Earthquake Sequence in Table II. The tabulated values correspond to AIC/Nc . The following conclusions can be reached from this Table: (i) The gamma, generalized Weibull, and κ-Weibull distributions have similar 20

TABLE I: p-values of Kolmogorov-Smirnov test for the fit between different probability models and the Crete Earthquake Sequence. Results are based on 1 000 simulations. The sample size (number of return intervals) is shown within parentheses next to the cutoff magnitudes. Gamma Weibull κ-Weibull Gen. Gamma Normal Lognormal Poisson ML,c =2.3 (628)

0.00

0.09

0.75

0.10

0.00

0.00

0.00

ML,c =2.5 (414)

0.07

0.40

0.14

0.35

0.00

0.00

0.00

ML,c =2.7 (273)

0.38

0.01

0.00

0.04

0.00

0.00

0.00

ML,c =2.9 (176)

0.14

0.23

0.20

0.25

0.00

0.00

0.00

ML,c =3.1 (103)

0.28

0.27

0.20

0.40

0.00

0.00

0.00

ML,c =3.3 (69)

0.28

0.11

0.10

0.19

0.00

0.01

0.00

ML,c =3.5 (45)

0.15

0.03

0.03

0.08

0.00

0.05

0.00

ML,c =3.7 (28)

0.16

0.10

0.08

0.07

0.00

0.00

0.14

ML,c =3.9 (18)

0.20

0.19

0.17

0.09

0.27

0.02

0.25

AIC values which are lower than the normal, lognormal, and Poisson models. (ii) The AIC values of the four top ranking distributions are quite close to each other. (iii) The κ-Weibull has the lowest AIC for the larger samples (i.e., those with Mc = 2.3, 2.5). The observation (i) also explains the somewhat unexpected outcome of Table I, namely, that the p-values of the generalized gamma and the κ-Weibull are not —for all magnitude cutoffs— equal or higher than the p-values of the respective subordinated distributions, i.e., the gamma and the Weibull respectively: The estimates of the probability model parameters are based on the minimization of the negative log-likelihood, which provides a different measure of the fit between the data and the model distribution than the Kolmogorov-Smirnov distance. We checked that the incongruence remains even if the likelihood optimization algorithm for the generalized gamma and the κ-Weibull is initialized by the respective optimal parameters of the gamma and Weibull distributions for the same data set. 21

TABLE II: Akaike Information Criterion (AIC) values per sample point for different probability models. The model parameters used maximize the likelihood of the Crete Earthquake Sequence given the respective model. The sample size (number of return intervals) is shown within parentheses next to the cutoff magnitudes. Gamma Weibull κ-Weibull Gen. Gamma Normal Lognormal Poisson ML,c =2.3 (628)

23.23

23.16

23.12

23.15

25.95

23.21

23.46

ML,c =2.5 (414)

24.05

24.02

24.00

24.02

26.40

24.12

24.29

ML,c =2.7 (273)

24.91

24.90

24.90

24.90

26.89

25.07

25.12

ML,c =2.9 (176)

25.79

25.78

25.79

25.79

27.59

25.91

26.00

ML,c =3.1 (103)

26.81

26.80

26.82

26.82

28.56

26.92

27.08

ML,c =3.3 (69)

27.58

27.59

27.62

27.62

29.23

27.76

27.86

ML,c =3.5 (45)

28.17

28.20

28.24

28.23

30.22

28.38

28.57

ML,c =3.7 (28)

29.30

29.36

29.43

29.41

30.67

29.79

29.46

ML,c =3.9 (18)

30.37

30.43

30.54

30.50

31.08

30.92

30.38

B.

Southern California Data

We also analyze an earthquake sequence which contains 2 446 events in Southern California (114◦ – 122◦ West longitude and 32◦ – 37◦ North latitude) down to depths of ≤ 20 (km) with magnitudes from 1 (except for 3 events at 0.5) up to 6.5 (ML ); 2 444 of these events have magnitudes less than 5.0 (ML ), whereas the two main shocks have 6.0 (ML ) and 6.5 (ML ). The events occurred during the period from January 1, 2000 until March 27, 2012 [73]. The p-values of the Kolmogorov-Smirnov test are listed in Table III. The Weibull and κ-Weibull distributions have practically the same p-values. The gamma and generalized gamma models, however, show overall better agreement with the observed return intervals than the Weibull or the κ-Weibull. Most of the p-values obtained for this data set are considerably lower than their counterparts for the Cretan data set. To ensure that this difference is not caused by an insufficient number of Monte Carlo simulations, we repeated the numerical experiment with 5 000 Monte Carlo simulations, which confirmed the results of Table III with minor changes in the p-values. On the other hand, as shown in Table IV, the 22

TABLE III: p-values of Kolmogorov-Smirnov test for the fit between different probability models and the Southern California sequence. Results are based on 5 000 simulations. The sample size (number of return intervals) is shown within parentheses next to the cutoff magnitudes. Gamma Weibull κ-Weibull Gen. Gamma Normal Lognormal Poisson ML,c =2.3 (1341) 0.0296

0.0000

0.0000

0.0000

0.0000

0.0000

0.0000

ML,c =2.5 (964)

0.0368

0.0000

0.0000

0.0000

0.0000

0.0000

0.0000

ML,c =2.7 (687)

0.0348

0.0000

0.0000

0.0000

0.0000

0.0000

0.0000

ML,c =2.9 (457)

0.0908

0.0002

0.0004

0.0000

0.0000

0.0000

0.0000

ML,c =3.1 (309)

0.5036

0.0638

0.0556

0.0440

0.0000

0.0000

0.0000

ML,c =3.3 (206)

0.0254

0.0158

0.0150

0.0270

0.0000

0.0000

0.0000

ML,c =3.5 (121)

0.1836

0.0110

0.0086

0.0152

0.0000

0.0000

0.0000

ML,c =3.7 (77)

0.0798

0.0304

0.0258

0.0412

0.0000

0.0002

0.0000

ML,c =3.9 (44)

0.0150

0.0450

0.0378

0.0580

0.0000

0.0964

0.0000

gamma, generalized gamma, Weibull and κ-Weibull distributions have similar AIC values. The low p-values are an indication that none of the models tested match the data very well in terms of the Kolomogorov-Smirnov distance. The gamma, generalized gamma, Weibull and κ-Weibull, however, are not rejected at the 1% level for most thresholds. It should be noted that recent arguments based on Bayesian analysis of hypothesis testing suggest that the significance level 0.05 used to reject the null hypothesis is overly conservative and should be shifter to 0.005 [74].

VI.

FIBER BUNDLE MODELS

Fiber Bundle Models (FBM) are simple statistical models that were introduced to study the fracture of fibrous materials [48]. To date they are used in many research fields, including fracture of composite materials [75], landslides [76], glacier avalanches [77] and earthquake dynamics [32, 43, 45]. In spite of their conceptual simplicity, FBMs exhibit surprisingly rich behavior. An FBM consists of an arrangement of parallel fibers subject to an external load F 23

TABLE IV: Akaike Information Criterion (AIC) values per sample point for the best-fit models of the Southern California sequence. The model parameters used maximize the likelihood of the Southern California sequence given the respective model. The sample size (number of return intervals) is shown within parentheses next to the cutoff magnitudes. Gamma Weibull κ-Weibull Gen. Gamma Normal Lognormal Poisson ML,c =2.3 (1341)

26.45

26.45

26.45

26.44

29.19

26.68

27.14

ML,c =2.5 (964)

26.98

27.01

27.01

27.01

29.73

27.26

27.79

ML,c =2.7 (687)

27.55

27.59

27.59

27.60

30.33

27.85

28.47

ML,c =2.9 (457)

28.36

28.41

28.42

28.40

30.92

28.65

29.29

ML,c =3.1 (309)

28.99

29.02

29.03

29.02

31.82

29.22

30.06

ML,c =3.3 (206)

29.28

29.30

29.31

29.31

32.69

29.43

30.88

ML,c =3.5 (121)

30.03

30.07

30.08

30.08

33.84

30.22

31.84

ML,c =3.7 (77)

30.64

30.72

30.81

30.72

34.65

30.89

32.74

ML,c =3.9 (44)

30.70

30.67

30.72

30.72

35.68

30.67

33.46

(Fig. 9). The fibers have random strength thresholds that represent the heterogeneity of the medium. Due to the applied loading, each fiber is deformed and subject to stress. If the stress applied to a specific fibre exceeds its failure threshold, the fiber ruptures and the excess load is redistributed either globally or locally between the remaining fibers. The ensuing redistribution of the load to the surviving fibers may trigger an avalanche of breaks. Each fibre break releases the elastic energy accumulated in the fibre.

A.

FBM Return Interval Statistics

We assume that the strain of the fiber bundle increases linearly with time t, i.e., ∝ t. Without loss of generality we set the elastic modulus, the initial length L and the strain rate equal to unity, and we use the elongation x instead of to measure the loading. The individual fibers have random failure thresholds xc with pdf fxc (x). Failed fibers are removed, and the stress is then redistributed between the surviving fibers using the equal load sharing rule. The energy of each avalanche is equal to the sum of the Hookean energies of the broken fibers [80]. Only events that exceed an energy threshold Ec are counted. The return intervals 24

F

L x

F FIG. 9: Schematic of a fiber bundle of initial length L elongated by x due to the loading force F . Broken fibers do not contribute to the strength of the bundle.

are measured as the time difference between two events with energy E ≥ Ec . Avalanches are considered to occur instantaneously. Fig. 10 illustrates the evolution of the avalanche events in time. The plots are obtained by loading a single bundle of 107 fibers the strength of which follows the Weibull pdf with m = 5, xs = 1. The avalanche sequences correspond to energy thresholds given by log10 Ec = 1, 2, where log10 is the logarithm with base 10. Figure 11 compares the quantiles of the return interval distribution obtained from a single bundle with those of the optimal κ-Weibull distribution for different energy thresholds. The optimal κ-Weibull parameters for each threshold are shown in Table V. The estimated κ values based on maximum likelihood are κ ˆ ≈ 2. This result suggests that κ > 1 values indicate a highly correlated system that can not be decomposed into RVEs. As stated in Subsection III A, the κ-Weibull pdf exhibits a power-law upper tail with exponent α = 1 + m/κ. For the FBM investigated above, the exponent of the return interval pdf is α ˆ ≈ [2.13, 2.20].

VII.

DISCUSSION AND CONCLUSIONS

We investigated the statistics of return intervals in systems that obey weakest-link scaling. We propose that the κ-Weibull distribution is suitable for finite-size systems (where the size is measured in terms of RVE size) and that the parameter κ is determined by the size of the 25

10

log (E)

6

log10Ec = 1

4 2 0.65

0.66

0.67

0.68

0.69

0.7

0.71

0.72

t

10

log (E)

6

log10Ec = 2

5 4 3 2 0.706 0.708 0.71 0.712 0.714 0.716 0.718 0.72 0.722 t

FIG. 10: Evolution of avalanche events in time for a bundle of 107 fibers with Weibull-distributed strength thresholds (m = 5, xs = 1). The most energetic avalanche (log10 Ef ≈ 7.31) is generated when the bundle fails (during the last event of the breaking sequence) at time tf ≈ 0.725.

system. A characteristic property of the κ-Weibull distribution is the transition from Weibull to power-law scaling in the upper tail of the pdf. This leads to an upper tail which decays slower than the tail of the respective Weibull distribution —a feature useful for describing the statistics of earthquake return intervals. The transition point depends on the system size and the Weibull modulus. Recent studies have identified a slope change in logarithmic plots of the ERI pdf, attributed to spatial (for earthquakes) or temporal (for lab fracturing experiments) nonstationarity of the background productivity rate [40]. We demonstrated that finite-size effects have a similar impact on the ERI pdf. Hence, finite size can explain deviations of earthquake return intervals from Weibull scaling without invoking non-stationarity (spatial or temporal) in the background earthquake productivity rate. 26

−5

10

log10Ec = 0.5

log10Ec = 1 −5

10 −6

10

−6

κ−Weibull

10 −6

−5

10

−6

10

−5

10

10

−3

−4

10

10

log10Ec = 1.5

log10Ec = 2

−4

10 −5

10

−5

10 −5

10

−4

−5

10

−4

10

10

−3

10

FBM

FIG. 11: Quantile-quantile plots of FBM return intervals (dimensionless) versus the best-fit κ-Weibull distribution. Data are based on a bundle of 5 × 107 fibers with strength thresholds drawn from the Weibull distribution with m = 5, xs = 1, unit elastic constant, and unit strain rate.

TABLE V: Maximum Likelihood estimates of κ-Weibull distribution parameters for the return intervals of Figure 11. log10 Ec

N

τˆs

m ˆ

κ ˆ

0.5

44094

1.2 × 10−6

2.4

2.1

1

8449

3.7 × 10−6

2.4

2.0

1.5

1686

1.0 × 10−5

2.6

2.2

2

311

3.1 × 10−5

2.6

2.3

27

In addition, we show that a distinct feature of the κ-Weibull distribution is the dependence of its hazard rate function: for m > 1 it increases with increasing time interval up to a certain threshold, followed by a ∝ 1/τ drop. This is in contrast with the Weibull hazard rate for m > 1 which increases indefinitely. Therefore, the κ-Weibull distribution allows for temporal clustering of earthquakes independently of the value of the Weibull modulus. The application of the κ-Weibull distribution to ERI assumes the following:

• Statistical stationarity, i.e., uniform ERI distribution parameters over the spatial and temporal observation window.

• Renormalizability of the interacting fault system into an ensemble comprising a finite number of independent effective RVEs with identical interval scale.

• Specific but simple functional form for the RVE survival probability given by (18).

We believe that the κ-Weibull distribution is also potentially useful for modeling the fracture strength of heterogeneous quasibrittle structures. The latter involve a finite number of RVEs and their fracture strength obeys weakest-link scaling [38, 46]. The connection between ERI power-law scaling and fracture mechanics pursued herein and in [11] also requires further research. Finally, we have used statistical methods (Kolmogorov-Smirnov test, Akaike Information Criterion) to compare different hypotheses for ERI distributions. For example, using the Kolmogorov-Smirnov test we showed that for one sequence of earthquakes the κ-Weibull, Weibull, gamma, and generalized gamma probability models are acceptable at the 5% level.

ACKNOWLEDGMENTS

M. P. Petrakis and D. T. Hristopulos acknowledge support by the project FIBREBREAK, Special Research Fund Account, Technical University of Crete. The Cretan seismic data were graciously provided by D. Becker, Institute of Geophysics, Hamburg University, Germany. We would also like to thank an anonymous referee for helpful suggestions. 28

Appendix A: Inverse Transform Sampling Method

We generate random numbers from the κ-Weibull and the generalized gamma distributions using the inverse transform sampling method. We first illustrate the algorithm for the κ-Weibull random numbers. d

1. We generate uniform random numbers uτ = U (0, 1). g(τ )

2. We employ the conservation of probability under the variable transformation τ → uτ , i.e., Fn (τ ) = FU (uτ ) = uτ ⇒ τ = Fn−1 (uτ ), n p m 2 2m 2 1 + a τ /n − aτ /n and a = τs−m . where Fn (τ ) =

(A-1)

3. The above in light of (A-1) leads to τ = Fn−1 (uτ ) =

n 1/m 2

τs u−1/n − u1/n τ τ

1/m

= τs [− lnκ (uτ )]1/m ,

κ = n−1 . (A-2)

In the case of the generalized gamma, the cumulative probability distribution is given by k m , (τ /τs ) , (A-3) Fn (τ ) = γ m where γ(α, x) is the incomplete gamma function defined by Z x 1 y α−1 γ(α, x) = dy α e−y/τs . Γ(α) 0 τs The random numbers τ are then given by inverting γ(α, x(τ )) = uτ , that is by 1/m τ = τs γn−1 (α; uτ ) .

[1] The term reverse Weibull refers to the distribution for maxima; this is known as the Weibull distribution in physical sciences and engineering. [2] Counterpart distributions with reversed supports are obtained for maxima. [3] E. J. Gumbel, Ann. Inst. Henri Poincar´e 5, 115 (1935). [4] W. Weibull, J. Appl. Mech. 18, 293 (1951). [5] D. Sornette, Critical Phenomena in Natural Sciences (Springer, Berlin, 2004).

29

[6] I. Eliazar and J. Klafter, Physica A 367, 106 (2006). [7] C. Franzke, Phys. Rev. E 85, 031134 (2012). [8] C. Nicolis and G. Nicolis, Phys. Rev. E 85, 056217 (2012). [9] H. L. D. de S. Cavalcante, M. Ori´ a, D. Sornette, E. Ott, and D. J. Gauthier, Phys. Rev. Lett. 111, 198701 (2013). [10] K. T. Tallakstad, R. Toussaint, S. Santucci, and K. J. M˚ aløy, Phys. Rev. Lett. 110, 145501 (2013). [11] D. T. Hristopulos and V. Mouslopoulou, Physica A 392, 485 (2013). [12] The terms interevent times, waiting times and recurrence intervals are also used. A subtle distinction can be made between recurrence intervals, which refer to events that take place on the same fault, and interocurrence intervals which encompass all faults in a specified region [32]. The statistical properties of recurrence intervals are more difficult to estimate, because less information is available for individual faults. The distinction is, nevertheless, conceptually important, since recurrence intervals characterize the one-body, i.e., the singlefault, problem, while interoccurrence intervals are associated with the activity of the manybody system [13]. Herein we use the term return intervals without further distinction. [13] D. Sornette, Phys. Rep. 313, 237 (1999). [14] J. Zhuang, D. Harte, M. J. Werner, S. Hainzl, source for Statistical Seismicity Analysis

and S. Zhou, Community Online Re-

(2012), 10.5078/corssa-79905851, available at

http://www.corssa.org. [15] P. Bak, K. Christensen, L. Danon, and T. Scanlon, Phys. Rev. Lett. 88, 178501 (2002). [16] A. Corral, Phys. Rev. E 68, 035102 (2003). [17] A. Saichev and D. Sornette, J. Geophys. Res. 112, B04313/1 (2007). [18] W. Klein, J. B. Rundle, and C. D. Ferguson, Phys. Rev. Lett. 78, 3793 (1997). [19] C. A. Serino, K. F. Tiampo, and W. Klein, Phys. Rev. Lett. 106, 108501 (2011). [20] A. Corral, Phys. Rev. Lett. 92, 108501 (2004). [21] A. Corral, Phys. Rev. Lett. 97, 178501 (2006). [22] A. Corral and K. Christensen, Phys. Rev. Lett. 96, 109801 (2006). [23] A. Saichev and D. Sornette, Phys. Rev. Lett. 97, 078501 (2006). [24] Y. Ogata, J. Am. Stat. Assoc. 83, 9 (1988). [25] Y. Hagiwara, Tectonophysics 23, 313 (1974).

30

[26] T. Rikitake, Tectonophysics 35, 335 (1976). [27] T. Rikitake, Tectonophysics 199, 121 (1991). [28] K. Sieh, M. Stuiver, and D. Brillinger, J. Geophys. Res. 94, 603 (1989). [29] G. Yakovlev, D. L. Turcotte, J. B. Rundle, and P. B. Rundle, Bull. Seismol. Soc. Am. 96, 1995 (2006). [30] J. R. Holliday, J. B. Rundle, D. L. Turcotte, W. Klein, K. F. Tiampo, and A. Donnellan, Phys. Rev. Lett. 97, 238501 (2006). [31] S. G. Abaimov, D. L. Turcotte, and J. B. Rundle, Geophys. J. Int. 170, 1289 (2007). [32] S. G. Abaimov, D. Turcotte, R. Shcherbakov, J. B. Rundle, G. Yakovlev, C. Goltz, and W. I. Newman, Pure Appl. Geophys. 165, 777 (2008). [33] T. Hasumi, T. Akimoto, and Y. Aizawa, Physica A 388, 483 (2009). [34] T. Hasumi, T. Akimoto, and Y. Aizawa, Physica A 388, 491 (2009). [35] M. S. Santhanam and H. Kantz, Phys. Rev. E 78, 051113 (2008). [36] D. T. Hristopulos and T. Uesaka, Phys. Rev. B 70, 064108 (2004). [37] S.-D. Pang, Z. Bazant, and J.-L. Le, Int. J. Fract. 154, 131 (2008). [38] Z. P. Bazant, J.-L. Le, and M. Z. Bazant, Proc. Natl. Acad. Sci. U.S.A. 1061, 11484 (2009). [39] P. M. Amaral, L. G. Rosa, and J. C. Fernandes, Rock Mech. Rock Eng. 41, 917928 (2008). [40] J. Bar´o, A. Corral, X. Illa, A. Planes, E. K. H. Salje, W. Schranz, D. E. Soto-Parra, and E. Vives, Phys. Rev. Lett. 110, 088702 (2013). [41] I. Main, Phys. 6, 20 (2013). [42] W. A. Curtin, Phys. Rev. Lett. 80, 1445 (1998). [43] M. J. Alava, P. K. V. V. Nukala, and S. Zapperi, Adv. Phys. 55, 349 (2006). [44] M. J. Alava, P. K. V. V. Nukala, and S. Zapperi, J. Phys. D 42, 214012 (2009). [45] B. K. Chakrabarti and L. G. Benguigui, Statistical Physics of Fracture and Breakdown in Disordered Systems (Clarendon Press, Oxford, UK, 1997). [46] Z. P. Bazant and S.-D. Pang, Proc. Natl. Acad. Sci. U.S.A. 103, 9434 (2006). [47] We assume a fixed RVE size with volume ∝ n. [48] H. E. Daniels, Proc. R. Soc. London Ser. A 183, 405 (1945). [49] R. L. Smith and S. L. Phoenix, J. Appl. Mech. 48, 75 (1981). [50] S. Phoenix, M. Ibnabdeljalil, and C.-Y. Hui, Int. J. Solids Struct. 34, 545 (1997). [51] I. Eliazar and R. Metzler, J. Chem. Phys. 137, 234106 (2012).

31

[52] I. Eliazar and R. Metzler, Phys. Rev. E 87, 022141 (2013). [53] Z. P. Bazant and S.-D. Pang, J. Mech. Phys. Solids 55, 91 (2007). [54] I. Eliazar and J. Klafter, Physica A 334, 1 (2004). [55] G. Kaniadakis, Eur. Phys. J. B 70, 3 (2009). [56] G. Kaniadakis, Phys. Lett. A 288, 283 (2001). [57] G. Kaniadakis, Phys. Rev. E 72, 036108 (2005). [58] G. Kaniadakis, Eur. Phys. J. A 40, 275 (2009). [59] G. Lapenta, S. Markidis, A. Marocchino, and G. Kaniadakis, Astrophys. J. 666, 949 (2007). [60] G. Kaniadakis, Entropy 15, 3983 (2013). [61] F. Clementi, M. Gallegati, and G. Kaniadakis, Eur. Phys. J. B 57, 187 (2007). [62] F. Clementi, T. Di Matteo, M. Gallegati, and G. Kaniadakis, Physica A 387, 3201 (2008). [63] F. Clementi, M. Gallegati, and G. Kaniadakis, J. Stat. Mech.: Theory Exp. P02037 (2009). [64] The number of RVEs n may also depend on the earthquake cutoff magnitude. [65] D. J. Daley and D. Vere-Jones, An Introduction to the Theory of Point Processes. Vol. I, 2nd ed. (Springer-Verlag, New York, 2003) pp. xxii+469. [66] D. Sornette and L. Knopoff, Bull. Seismol. Soc. Am. 87, 789 (1997). [67] A. Corral, Phys. Rev. E 71, 017101 (2005). [68] S. Touati, M. Naylor, and I. G. Main, Phys. Rev. Lett. 102, 168501 (2009). [69] We use the MATLAB fmincon constrained minimization function with a trust region reflective algorithm and explicit gradient information to minimize the negative log-likelihood. The method is applied iteratively (i.e., fmincon is called with the last estimate as initial guess) for fifty times. The optimization parameters used are: Maximum number of objective function evaluation = 2 × 104 , maximum number of iterations = 2 × 104 , objective function error tolerance = 1 × 10−5 . Tests with the data sequences showed that this number of iterations is sufficient for the parameter estimates to converge. The accuracy and precision of the estimates were also tested using synthetic sequences of κ-Weibull random numbers generated by means of the inversion method. [70] A. Clauset, C. Shalizi, and M. Newman, SIAM Rev. 51, 661 (2009). [71] E. W. Stacy, Ann. Math. Stat. 33, 1187 (1962). [72] H. Akaike, IEEE Trans. Automat. Contr. 19, 716 (1974).

32

[73] The facilities of the Southern Californian earthquake Data Center (SCEDC), and the Southern California Seismic Network (SCSN), were used for access to waveforms, parametric data, and metadata required in this study. The SCEDC and SCSN are funded through U.S. Geological Survey Grant G10AP00091, and the Southern Californian earthquake Center, which is funded by NSF Cooperative Agreement EAR-0529922 and USGS Cooperative Agreement 07HQAG0008. [74] V. E. Johnson, Proc. Nat. Acad. Sci. 110, 19313 (2013). [75] S. L. Phoenix and L.-J. Tierney, Eng. Fract. Mech. 18, 193 (1983). [76] D. Cohen, P. Lehmann, and D. Or, Water Resour. Res. 45, W10436 (2009). [77] I. Reiweger, J. Schweizer, J. Dual, and H. J. Herrmann, J. Glaciol. 55, 997 (2009). [78] B. Gutenberg and C. F. Richter, Bull. Seismol. Soc. Am. 46, 105 (1956). [79] B. Gutenberg and C. F. Richter, Seismicity of the Earth and Associated Phenomena (Hafner New York, 1965). [80] S. Pradhan and P. C. Hemmer, Phys. Rev. E 77, 031138 (2008).

33