On the earthquake predictability of fault interaction models.

Space-time clustering is the most striking departure of large earthquakes occurrence process from randomness. These clusters are usually described ex-post by a physics-based model in which earthquakes are triggered by Coulomb stress changes induced by other surrounding earthquakes. Notwithstanding the popularity of this kind of modeling, its ex-ante skill in terms of earthquake predictability gain is still unknown. Here we show that even in synthetic systems that are rooted on the physics of fault interaction using the Coulomb stress changes, such a kind of modeling often does not increase significantly earthquake predictability. Earthquake predictability of a fault may increase only when the Coulomb stress change induced by a nearby earthquake is much larger than the stress changes caused by earthquakes on other faults and by the intrinsic variability of the earthquake occurrence process.

1. Introduction

Since early 1990s, the space-time earthquake clustering has been mostly explained by fault interaction [Stein et al., 1992; King et al., 1994]. When an earthquake occurs, the state of stress in the Earth is altered as a consequence of the seismic dislocation. The change in the stress field can be predicted by means of standard elasticity theory if the parameters describing the seismogenic fault are known, and its effect on nearby faults can be quantitatively assessed through the evaluation of the so-called Coulomb Failure Function (CFF) [Reasenberg and Simpson, 1992]. If we indicate with Δτ and Δθ the changes in shear and normal stresses relative to a given fault plane and with Δp the pore pressure change, the variation of the CFF is defined as ΔσCFF=Δτ + μ(Δp−Δθ), where μis the friction coefficient on the fault. A rupture is promoted by increasing the shear (Δτ > 0) and/or fault unclamping (Δθ < 0), while the effect of a pore pressure increase to counteract normal stress is often accounted by defining an effective friction coefficient μ′ and writing ΔσCFF=Δτ − μ′Δθ. The effective friction coefficient depends on static friction, hydrostatic pressure, and pore fluid pressure; a value of μ′ = 0.4, consistent with laboratory evidences on friction and moderate pore pressure when fluids are not fully expelled, is often assumed [Stein et al., 1992, 1997]. This model has been later generalized to account for the additional viscoelastic response of the Earth's layers [Pollitz, 1992; Piersanti et al., 1997]. In essence, this element brings to an extra contribution of ΔσCFF that is accumulated through time.

The modeling based on CFF (CFF-modeling hereafter) became immediately very popular, and it got a wide acceptance from the scientific community, as testified by the plethora of papers that report retrospective analysis of large earthquakes cluster in terms of CFF-modeling [e.g., Stein, 1999; Steacy et al., 2005; Nalbant et al., 2005; Pollitz et al., 2006; Stramondo et al., 2011]. However, the very few efforts devoted to use CFF-models to forecast earthquakes in a prospective way did not lead to particularly encouraging results [Woessner et al., 2011; Parsons et al., 2012; Steacy et al., 2014]. It has been argued that the negligible gain in earthquake predictability of CFF-models may be due to the limited knowledge of the faults geometry and other source parameters that can prevent the real-time application of such a kind of models [Steacy et al., 2014]. For this reason, in this paper we explore the earthquake predictability of CFF-models in an ideal case, where all relevant parameters are perfectly known in advance.

2. The Model and the Synthetic Catalogs

We analyze synthetic catalogs generated by a model entirely based on CFF-modeling. The model used here is purposely simple in order to better control the effects of each model component on the results. The model has three main ingredients. First, we consider a realistic set of faults; specifically, we use the real faults distribution of Central Italy (Figure 1, fault network C-ITALY hereafter). The faults, their geometry and parameters, are taken from Marzocchi et al. [2009, and references therein], and we assume that they are all perfectly known. Second, we impose that each fault behaves like a Brownian Random Oscillator (BRO) [Matthews et al., 2002] and produces only one characteristic magnitude [Schwartz and Coppersmith, 1984]. In this model, the stress in each fault is loaded as a Brownian motion and it is described by a drifted Wiener process [Matthews et al., 2002], where a constant tectonic drift is perturbed by random fluctuations; these fluctuations mimic the influence of several factors, like the stress perturbations induced by external sources, the stress weakening or hardening of the fault, pore pressure variations, etc. Once the stress achieves a critical value δ, an earthquake occurs and sets the stress to zero; the fluctuations of the stress loading may also accommodate a possible variability of δ. In our application, we set the failure threshold δ = 3 MPa that is a reasonable average value for normal faults in Italy [Cocco and Rovelli, 1989] and at global scale [Allmann and Shearer, 2009]. Matthews et al. [2002] show that this conceptual BRO model leads to a Brownian Passage Time (BPT) distribution for the interevent times between consecutive earthquakes. This distribution is governed by two parameters, the average interevent times and the coefficient of variation α(also called aperiodicity) that is the ratio between the standard deviation and the mean of the interevent times. Here we set years, and α = 0.3. Usually, the aperiodicity observed on faults is larger than this value. However, to avoid a double counting of fault interaction effects, we use the aperiodicity that has been observed on one isolated fault, i.e., a fault that does not have significant interactions with other faults [Berryman et al., 2012]; in this case, the aperiodicity αdescribes only the intrinsic earthquake occurrence randomness. The third ingredient of the model is given by fault interaction, i.e., each earthquake produces a stress change in the nearby faults that follows the CFF principles. The coseismic and postseismic stress changes are modeled through equations (3) and (4) in Marzocchi et al. [2009]. We assume a uniform slip on the source fault, and the stress change on the receiver fault is the average stress calculated over the whole fault plane. From a physical point if view, the use of the average stress change implies that the nucleation point in the receiver fault has a uniform spatial distribution, and it changes location from earthquake to earthquake.

Figure 1

Map of the major seismogenic faults in Central Italy [Basili et al., 2008; Marzocchi et al., 2009], named C-ITALY fault network. The red boxes represent the set of faults that mimics a strongly coupled simplified network geometry (LINE fault network). The two target faults are marked by the abbreviations OPF and FF.

Once the model has been defined, we generate a synthetic seismic catalog 108 years long in order to get a reasonable statistics on earthquake occurrences. In particular, we focus our attention on the analysis of earthquake occurrences associated with two target seismogenic structures, Fucino and Ovindoli Pezza faults (hereafter F-fault and OP-fault, respectively), that are among the largest known in Central Italy (Figure 1). These two faults have a strong positive coupling in terms of CFF (Figure 2). An earthquake on F-fault (M = 6.7) induces a ΔσCFF∼1.6 × 10−1MPa on OP-fault (Figure A1 in the supporting information); this value is larger than the range of 10−1−10−2 MPa that characterizes the typical stress transfer between faults [Scholz, 2010]. Here we aim at exploring whether fault interaction modifies the earthquake predictability skill on OP-fault and/or brings to a time synchronization of the two target faults.

Figure 2

Empirical cumulative distribution of the stress changes induced by the C-ITALY fault network on OP-fault. The red dot is the F-fault contribution.

3. Earthquake Predictability and Synchronization of Synthetic Catalogs

We investigate on earthquake predictability of CFF-modeling comparing two different probabilities for OP-fault. First, we consider the probability, Pr(EOPF|ΔtOPF), of an earthquake on OP-fault given the elapsed time ΔtOPF since the previous earthquake on the same fault. This probability is evaluated numerically by the frequency n/N, where n is the number of earthquakes on the OP-fault with elapsed time in the time bin [ΔtOPF−25,ΔtOPF+25] years, and N is the total number of earthquakes on the OP-fault regardless of the elapsed time. The amplitude of the time bin is chosen to account for two opposite requirements: it should be long enough to filter short-term fluctuation of the process (two elapsed times are considered identical if their difference is smaller than ΔtOPF), but not too long to blur the time dependency of the process, i.e., . Second, we consider the probability Pr(EOPF|ΔtOPF,EFF) as the frequency n′/N′, where N′ is composed by a subset of N containing only earthquakes on the OP-fault that were anticipated by earthquakes on the F-fault in the last 10 years. This time interval describes well the time clustering between large earthquakes [Parsons, 2002; Faenza et al., 2008]. As before, the variable n′ is the number of earthquakes in N′ with elapsed time in the time bin [ΔtOPF−25,ΔtOPF+25] years.

The ratio of these two probabilities, G = Pr(EOPF|ΔtOPF,EFF)/Pr(EOPF|ΔtOPF), represents the probability gain due to the knowledge of fault interaction. We calculate the uncertainty associated to this ratio, assuming that the variables n, n′, and N′ are counts from the Poisson process, where the standard deviation is the square root of the count. Then, the error bar of the two frequencies and of the ratio G can be obtained propagating these standard deviations. Figure 3a shows the ratio of these two probabilities G, and its standard deviation, σ. Note that G is never significantly larger than 1. This result implies that even in the ideal case in which CFF represents the real physics of fault interaction, the knowledge of an earthquake that occurred on F-fault does not improve significantly the forecasting capability on OP-fault that can be obtained by the earthquake statistics observed on the OP-fault.

Figure 3

Results of the model using C-ITALY and α = 0.3. (a) Probability gain G = Pr(EOPF|ΔtOPF,EFF)/Pr(EFF|ΔtOPF); the horizontal dotted line marks the case G = 1 of no probability gain. The vertical bar shows the one-sigma uncertainty σ. (b) Probability density function (PDF) of the interevent times between events on OP-fault and F-fault; the red line shows the distribution when the earthquakes occur randomly on OP-fault and F-fault (f(x(rand)) in the text), and the black line shows the distribution observed in the synthetic catalogs (f(x(obs)) in the text). The vertical dashed lines show the average of these distributions.

Fault interaction can bring to a time synchronization of faults [Scholz, 2010]. In this case, the probability gain due to the knowledge of an earthquake on F-fault may be negligible, simply because the synchronization implies that an earthquake on F-fault occurs when OP-fault is very close to the end of its natural cycle; this means that these two information (the occurrence of an earthquake on F-fault and the time since the last earthquake on OP-fault) are redundant. We investigate the time synchronization of faults analyzing the statistical distribution of interevent times, i.e., , where is the time occurrence of the ith earthquake on the OP-fault, is the time occurrence of the closest in time earthquake on the F-fault, and . In particular, we calculate the difference where E[x(obs)] is the expectation value of the distribution f(x(obs)) of the interevent times that are observed in the synthetic catalog and E[x(rand)] is the expectation of f(x(rand)) that is calculated when the earthquakes occur randomly on OP-fault and F-fault. D close to 1 means that there is a marked decrease of the interevent times (E[x(rand)] >> E[x(obs)]) that is typical for a time synchronization of faults, while a small D means that fault interaction does not modify the central value of f(x(rand)). Figure 3b shows that the differences between the two distributions and their central values are negligible.

The negligible gain in earthquake predictability (Figure 3a) and in time synchronization (Figure 3b) can be explained by analyzing the relationship between the ΔσCFF induced by F-fault on OP-fault and the intrinsic variability of the stressing rate. The BPT distribution is based on the modeling of the stressing rate by means of a Wiener process [Matthews et al., 2002], in which the stress changes (ΔσCFF) have amplitude that is related to the coefficient of variation αof the BPT distribution; in particular, the variance of the stress changes of the Wiener process is given by var(ΔσBPT) = (αδ)2[Matthews et al., 2002]. In our model, to this intrinsic variability of the BRO model, we have to sum the variability of the independent stress contributions given by all other faults (ΔσCFF); so assuming that the variability of the BRO model and ΔσCFF induced by the other faults are independent, the standard deviation of the overall stress changes is given by where var(ΔσCFF) has been evaluated numerically by the simulations as the variance of the ΔσCFF induced by all faults on F-fault in a time interval of . We evaluate quantitatively the effect of F-fault on OP-fault using the influence parameter ξ that is defined as where the asterisk marks the calculation made without the F-fault contribution. The influence parameter ξ is the normalized difference of the standard deviations when the effect of the F-fault on OP-fault is included and when it is not. In our model ξ = 3×10−4.

The results of Figure 3 show that even in a simple and perfectly known physical system, the earthquake predictability of a fault, based on empirical statistics of past earthquakes, is not increased by the knowledge of the Coulomb stress change derived from nearby earthquakes. In order to explore the conditions that may lead to some increase of earthquake predictability, we rerun the same calculations using different parametrizations of the model. In particular, we analyze the importance of the fault network complexity by considering a simplified fault network that consists only of few aligned and strongly coupled faults (the red faults in Figure 1, fault network LINE hereafter). We evaluate the role of the intrinsic randomness of the earthquake occurrence process using δ = 2 MPa, α = 0.15, and α = 0. All these parameters lead to a more deterministic process; in particular, α = 0means that the BRO model becomes a pure deterministic oscillator and ψ is given only by the standard deviation of the ΔσCFF induced by all faults on OP-fault.

Figure 4 and Table 1 show the results for different model parametrizations (see also Figures A2–A5 in the supporting information). Probability gain is almost never significantly larger than 1; there is only one parametrization (with α = 0, C-ITALY fault network) in which the probability gain is significantly larger than 1, and it is about 5. The synchronization parameter shows that for most parametrizations, D is small and faults do not synchronize in time. Just in one parametrization (LINE fault network and α = 0), D is close to 1, implying strong time synchronization. All these behaviors are well described by the influence factor ξ. For very small values of ξ, the probability gain and the time synchronization are both negligible. The largest ξ values (∼0.3) are obtained for strongly deterministic processes (α = 0) where OP-fault and F-fault have a strong time synchronization (D∼ 1) or show a marked probability gain. In summary, the joint effect in fault interaction modeling played by the fault network geometry, the intrinsic variability of the earthquake occurrence process on each fault, and magnitude of the stress changes induced by surrounding earthquakes can be captured by the influence parameter that is positively correlated to the increase of earthquake predictability of a fault.

Figure 4

As for Figure 3 but relative to some different parametrizations of the model. The results for other parametrizations of the model are reported in the supporting information.

Table 1

Results for Different Parametrizations of the Modela

Parametrization of the Model z	Probability Gain at Maximum	Synchronization Parameter D	Influence Factor ξ
C-ITALY, α = 0.3, δ = 2MPa	1.3 ± 0.2	5.2 × 10−2	2.8 × 10−4
C-ITALY, α = 0.3, δ = 3MPa	1.3 ± 0.2	6.4 × 10−2	3.0 × 10−4
C-ITALY, α = 0.15, δ = 3MPa	1.3 ± 0.2	2.8 × 10−3	1.2 × 10−3
C-ITALY, α = 0.15, δ = 2MPa	1.1 ± 0.1	6.3 × 10−3	1.2 × 10−3
LINE, α = 0.3, δ = 2MPa	1.2 ± 0.1	4.6 × 10−2	1.3 × 10−3
LINE, α = 0.3, δ = 3MPa	1.4 ± 0.2	6.1 × 10−2	1.3 × 10−3
LINE, α = 0.15, δ = 2MPa	1.1 ± 0.1	1.7 × 10−3	5.2 × 10−3
LINE, α = 0.15, δ = 3MPa	1.1 ± 0.1	1.4 × 10−3	5.2 × 10−3
LINE, α = 0, δ = 3MPa	1.0 ± 0.1	9.9 × 10−1	2.9 × 10−1
C-ITALY, α = 0, δ = 3MPa	4.9 ± 0.3	2.3 × 10−4	3.0 × 10−1

The second column shows the probability gain G having the maximum normalized value z = (G − 1)/σ. The cases are ordered for increasing values of the influence factor ξ.

Finally, we have checked the stability of the results including the postseismic effect as modeled by Marzocchi et al. [2009] and using a much shorter interevent times average ( years) for all faults. The results remain essentially the same (see Figures A4 and A5 in the supporting information); the largest discrepancy is discussed in the caption of Figure A4.

4. Discussion and Conclusions

The physical system adopted in this paper is a simplified version of the reality, but it has been designed to account for in a proper way some important aspects of the real system, such as the geometry of the fault network [Marzocchi et al., 2009], the intrinsic variability of the earthquake occurrence process [Matthews et al., 2002; Berryman et al., 2012], and the CFF-modeling [e.g., Stein, 1999]. We find that even in this simple system, the influence parameter is small and the increase in earthquake predictability using CFF-modeling is negligible. This is due to the fact that the CFF stress changes induced by F-fault on OP-fault have a minor effect with respect to the intrinsic randomness of the earthquake generation process; here the state of stress on OP-fault is continuously shuffled, and it is kept far from significant time synchronization with any adjacent fault. Worthy of note, this marked randomness remains even for a more deterministic configuration of the model with α = 0.15 and δ = 2 MPa (Table 1 and Figure 4).

We find that the CFF-modeling can increase the earthquake predictability on one fault only when the CFF stress change from a specific source fault is largely predominant over the stress changes caused by the other remaining faults and over the intrinsic variability of the earthquake occurrence process. In this case, we may either observe a marked time synchronization between faults or a significant probability gain; i.e., the occurrence of an earthquake on F-fault increases significantly the probability to have a close in time earthquake on OP-fault.

These results have an important physical consequence; even though CFF-modeling is able to explain the space-time earthquake clusters, it cannot help to foresee which fault will be the next one to slip. We obtain these results for an optimal case, where the physical model is perfectly known. The reality is certainly more complex, and our knowledge is limited; so it is expected that the real forecasting skill of CFF-modeling and the time synchronization of interacting faults are even less than what shown here.