Order parameter fluctuations in natural time and b-value variation before large earthquakes

Self-similarity may stem from two origins: the process' increments infinite variance and/or process' memory. The $b$-value of the Gutenberg-Richter law comes from the first origin. In the frame of natural time analysis of earthquake data, a fall of the b-value observed before large earthquakes reflects an increase of the order parameter fluctuations upon approaching the critical point (mainshock). The increase of these fluctuations, however, is also influenced from the second origin of self-similarity, i.e., temporal correlations between earthquake magnitudes. This is supported by observations and simulations of an earthquake model.

A large variety of natural systems exhibit irregular and complex behavior which at first look seems to be erratic, but in fact possesses scale-invariant structure, for example see Refs. [1,2]. A stochastic process X (t) is called self-similar [3] with index H > 0 if it has the property where the equality concerns the finite-dimensional distributions of the process X (t) on the rightand the left-hand side of the equation (not the values of the process).
A point of crucial importance in analyzing data from complex systems that exhibit scaleinvariant structure, is the following: In several systems this nontrivial structure stems from longrange temporal correlations; in other words, the self-similarity originates from the process' memory only. This is the case for example of fractional Brownian motion. Alternatively, the selfsimilarity may solely come from the process' increments infinite variance. Such an example is Lévy stable motion (the variance of Lévy stable distributions is infinite since they have heavy tails [4], thus differing greatly from the Gaussian ones). In general, however, the self-similarity may result from both these origins [5], the presence of which can be in principle identified when analyzing the complex time series in terms of the new time domain termed natural time [6].
The evolution of seismicity is a typical example of complex time series. Several traditional studies were focused on the variation of the b-value of the Gutenberg-Richter (G-R) law [7], which states that the (cumulative) number of earthquakes with magnitude greater than (or equal to) M, N(≥ M), occurring in a specified area and time is given by where b is a constant, varying only slightly from region to region and the constant a gives the logarithm of the number of earthquakes with magnitude greater than zero [8]. These studies found that the b-value decreases before a large event, e.g., see Ref. [9] (cases where b-value increases prior to and then decreases sharply before a large event have been also reported [10]). Here, considering that the b-value itself solely focuses on the one origin of self-similarity, and in particular the process' increments infinite variance, we show that, when employing natural time analysis, the b-value decrease before large earthquakes reflects an increase of the fluctuations of the order parameter of seismicity when approaching the critical point (mainshock, see below). The whole precursory variation of the order parameter fluctuations, however, is more complex since it captures both origins. Temporal correlations between earthquake magnitudes also play an important role in this precursory variation, thus leading to more spectacular results compared to the ones obtained when restricting ourselves to traditional analysis of b-value alone.
For a time series comprising N events, we define [11] the natural time χ k for the occurrence of the k-th event (of energy Q k ) by χ k = k/N. We then study the evolution of the pair (χ k , Q k ) or where ω stands for the natural angular frequency, and then evaluate the real function Π(ω) = |Φ(ω)| 2 in the low frequency limit. By considering the Taylor expansion Π(ω) = 1 − κ 1 ω 2 + κ 2 ω 4 + . . . , we find that the approach of a dynamical system to criticality (see Chapter 8 of Ref. [6]) is identified by means of κ 1 , i.e., which is the variance [6,11,12] of natural time weighted for p k . When Q k are independent and identically distributed positive random variables, we obtain the "uniform" (u) distribution of p k , as it was defined in Ref. [13] (see also p.122 of Ref. [6]). In this case, all p k vary around their mean value 1/N (cf. since ∑ N n=1 p n = 1) and the quantity κ 1 results [13] in κ u = 1/12 for large N. In general, in a complex time series, in order to identify the two origins of self-similarity by means of natural time analysis, we focus on the expectation value E (κ 1 ) of the variance κ 1 of natural time when sliding a natural time window of length l through a time series of Q k > 0, k = 1, 2, . . .N.
If self-similarity exclusively results from the process' memory, the E (κ 1 ) value should change to κ u = 1/12 for the (randomly) shuffled data. This is the case of the Seismic Electric Signals (SES) activities [14], which are series of low-frequency (≤ 1Hz) electric signals detected a few to several weeks (up to five months) before an earthquake when the stress in the focal region reaches a critical value (and hence long range correlations develop). For example, the three upper channels in Fig.1(b) show three SES activities that preceded major earthquakes in southern, southwestern and western Greece, respectively, as depicted in the map of Fig.1(a). For the sake of comparison, the lowest channel shows an SES activity recorded in northern Greece (close to Thessaloniki). In all these four cases, the analysis of their original data lead to κ 1 ≈ 0.07 (see also below), which turns to κ u = 1/12 upon shuffling the data. On the other hand, if the self-similarity results from process' increments infinite variance only, E (κ 1 ) should be the same (but differing from κ u ) for the original and the (randomly) shuffled data. Finally, when both origins of self-similarity are present, the relative strength of the contribution of the one origin compared to that of the other can be quantified on the basis of Eqs. (12) and (13) of Ref. [15] (see also Ref. [6]).
In what remains, we focus on complex time series of seismicity. Earthquakes exhibit scaling relations chief among which is the aforementioned G-R law [7]. For reasons of convenience, we write hereafter G-R law of Eq.(2) into the form N(≥ M) ∝ 10 −bM . Considering that the seismic energy E released during an earthquake is related [16] to the magnitude through E ∝ 10 cM , where c is around 1.5, the latter form turns to the distribution, where γ = 1 + b/1.5. Hence, b ≈ 1 means that the exponent γ is around γ=1.6 to 1.7, see Table 2.1 of Ref. [6].
The complex correlations in time, space and magnitude of earthquakes have been extensively studied [17][18][19][20][21]. The observed earthquake scaling laws [22] seem to indicate the existence of phenomena closely associated with the proximity of the system to a critical point (e.g., see Ref. [18] and references therein). In the frame of natural time analysis, it has been suggested [12] (see also pp.249-254 of Ref. [6]) that the order parameter of seismicity is the quantity κ 1 . The κ 1 value itself may lead to the determination of the occurrence time of the impending mainshock [6,11,15,23] when SES data are available. In particular, when the κ 1 value resulting from the natural time analysis of the seismicity subsequent to the SES recording becomes approximately equal to 0.070, the mainshock occurs within a time window of the order of one week. This has been empirically observed in several cases [11,15,23] (see also Chapter 7 of Ref. [6]) including the three major earthquakes of Fig.1(a) that followed the SES activities depicted in Fig.1(b). An example of the κ 1 dynamics after the recording of the SES activity depicted in the third channel of Fig.1(b) until the occurrence of the magnitude 6.4 mainshock on June 8, 2008 (blue star in Fig.1(a)) is given in Ref. [24]. In the lack of SES data, we have to solely rely on the fluctuations of the order parameter of seismicity. Along these lines, we investigated [25] the period before and after a significant mainshock. Time-series for various lengths of W earthquakes that occurred before or after the mainshock have been studied. The probability distribution function (pdf) P(κ 1 ) versus κ 1 was found to exhibit a bimodal feature when approaching a mainshock. To quantify this feature, we considered the variability of κ 1 , which is just the ratio where σ (κ 1 ) and µ(κ 1 ) stand for the standard deviation and the mean value of κ 1 for sliding window lengths l=6-40. The bimodal feature reflects that, upon approaching the mainshock (with the number W of the earthquakes before mainshock decreasing), the variability of κ 1 should increase. This was subsequently confirmed because before the M9.0 devastating Tohoku earthquake in Japan on March 11, 2011, the variability of κ 1 exhibited [26] a dramatic increase.
In addition, we investigated [27] Fig.2(b). The general feature of this curve is more or less similar to that observed for example before Tohoku earthquake [26]; quantitative agreement cannot be demanded, however, because temporal correlations between the earthquake magnitudes are also present which influence the observed results. This is corroborated by the following results obtained from the Olami-Feder-Christensen (OFC) earthquake model [28]. We preferred to employ this model here, since it has been studied in detail in hundreds of publications, but we clarify that there exist more recent ones, e.g., see Ref. [29] where the primary role of the fault system geometry is emerged.
The OFC model runs as follows: we assign a continuous random variable z i j ∈ (0, 1) to each site of a square lattice, which represents the local "energy". Starting with a random initial configuration taken from a uniform distribution in the segment (0,1), the value z i j of all sites is simultaneously increased at a uniform loading rate until a site i j reaches the threshold value z thres =1 (i.e., the loading ∆ f is such that z i j max + ∆ f = 1). This site then topples which means that z i j is reset to zero and an "energy" αz i j is passed to every nearest neighbor, where the coupling parameter α can take values from zero to 0.25 and is the only parameter of the model, apart from the edge length L of the square lattice. If this causes a neighbor to exceed the threshold, the neighbor topples also, and the avalanche continues until all z kl < 1. Then the uniform loading increase resumes. The number of topplings defines the size of an avalanche or "earthquake" and (when it is larger than unity k increases by one) is used as Q k in natural time analysis. Here, we use the case of free boundary conditions [30] in which α varies locally α i j = 1 n i j +K , where n i j is the actual number of nearest neighbors of the site i j (for sites in the bulk n i j = 4, for sites at the edges n i j = 3 and for the four sites at the corners n i j = 2) and K denotes [30] the elastic constant of the upper leaf springs measured relatively to that of the other springs between blocks in the Burridge-Knopoff model [31]. The OFC model is obviously non-conservative for K > 0 for which α i j < 0.25 in the bulk (for more details on the OFC modelling see pp. 349-363 of Ref. [6] and references therein).
We first study the predictability of the OFC model on the basis of the κ 1 variability. We consider the variability β k which is a function of the natural time index k, k = 1, 2, . . ., N = 2 × 10 6 estimated by analyzing in natural time for each k the preceding W =100 avalanches. The time increased probability (TIP) [32] (i.e., the time during which there exists a high probability for the occurrence of a large avalanche exceeding a given threshold) is turned on when β k > β c , where β c is a given threshold in the prediction. If the size Q k is greater than a target avalanche size threshold Q c , we have a successful prediction. For binary predictions, the prediction of events becomes a classification task with two types of errors: missing an event and giving a false alarm. We therefore choose [33] the receiver operating characteristics (ROC) graph [34] to depict the prediction quality. This is a plot of the hit rate versus the false alarm rate, as a function of the total rate of alarms, which here is tuned by the threshold β c . Only if in between the hit rate exceeds the false alarm rate, the predictor is useful. Random predictions generate equal hit and alarm rate, and hence they lead to the diagonal in ROC plot. Thus, only when the points lie above this diagonal the predictor is useful. As an example, the ROC graphs for L = 512 and K = 1 or L = 256 and K = 2 are shown in Fig. 3 (the rational for choosing these two cases stems from the study of Ref. [35] in which it was shown that the OFC model with free boundary conditions exhibits in these cases -see their shuffled and then β k was estimated (magenta curves); in both cases, we obtain curves which almost coincide with the diagonal. This clearly demonstrates that the aforementioned excess of the results related with the original Q k series from the diagonal comes from the sequential order of avalanches and cannot be considered as chancy.
We now proceed to the investigation of the temporal correlations between the magnitudes m k = log 10 (Q k )/1.5 obtained from the sizes Q k of the avalanches in the OFC model preceding a large avalanche. The results can be visualized in two examples in Fig.4, where we plot in blue the exponent a DFA of the detrended fluctuation analysis (DFA) [36] (along with the variability β plotted in red) versus the number W of avalanches before a large avalanche (negative x semi-axis, x = −W ). Note that DFA has already been employed in Ref. [37] for monitoring temporal correlations before bifurcations. In the upper example, Fig.4(a), the value of a DFA well before the large avalanche, being somewhat larger than 0.5, exhibits small changes but strongly increases upon approaching the large avalanche, i.e., at W = 100 the value of a DFA becomes ≈ 0.75 which shows intensified temporal correlations. In the lower example, Fig.4(b), well before the large avalanche we have a DFA ≈ 0.6 showing long range temporal correlations, which first turn to anticorrelations upon approaching the large avalanche, e.g., a DFA ≈ 0.43 at W = 400, and finally become random, i.e, a DFA ≈ 0.5 at W = 100, just before the "mainshock". In both examples of large (Q k > 30, 000) avalanches, only in 30% of the cases a rapid increase of β upon approaching them is observed. This is more or less consistent with empirical observations since in Japan this precursory increase was observed in 8 out of 25 earthquakes (all above M7 during 1 January 1994 to 11 March 2011 with depths smaller than 700 km) [26]. Concerning the α values, when studying W = 100, 200, . . .1000, among the 579 large avalanches studied, in 76% of the cases the α value was found to become smaller than 0.5 (as seen in Fig.4(b)). * Correspondence to: P. Varotsos (pvaro@otenet.gr) [1] C. K. Peng