- Email: [email protected]
Non-Critical Precursory Accelerating Seismicity Theory (NC PAST) and limits of the power-law fit methodology
Tectonophysics 452 (2008) 42–50
Contents lists available at ScienceDirect
Tectonophysics j o u r n a l h o m e p a g e : w w w. e l s ev i e r. c o m / l o c a t e / t e c t o
Non-Critical Precursory Accelerating Seismicity Theory (NC PAST) and limits of the power-law fit methodology A. Mignan Risk Management Solutions, Science and Technology Research, Peninsular House, 30 Monument Street, London EC3R 8NB, United Kingdom
a r t i c l e
i n f o
Article history: Received 18 September 2007 Received in revised form 2 February 2008 Accepted 13 February 2008 Available online 4 March 2008 Keywords: Earthquake Forecast Precursory seismicity AMR Stress loading Stress shadow Synthetic catalogue
a b s t r a c t The hypothesis that Accelerating Moment Release (AMR) is a precursor to large earthquakes is still debated. On one hand, AMR has been claimed to be observed in many cases and on the other hand, it has been proposed that apparent AMR is only due to data-fitting. The debate is in general focused on the validity of the c-value (curvature parameter), which permits to quantify AMR (i.e. cumulative Benioff strain through time), or more generally Precursory Accelerating Seismicity (PAS, i.e. cumulative number of events through time). Contrary to previous studies, which compare c-value optimization in real seismicity catalogues and in random synthetic catalogues, I test c-value optimization in theoretical synthetic catalogues. In that particular case, I assume that PAS exists and that it can be explained by the Non-Critical Precursory Accelerating Seismicity Theory (NC PAST). This theory demonstrates that PAS can emerge from the background seismicity because of the decrease, due to loading, of the size of a stress shadow due to a previous earthquake. I improve the NC PAST by integrating the following characteristics of the background seismicity, (1) the density of random events outside the stress shadow δ0b and (2) the noise ratio δ−b/δ0b, with δ−b being the density of random events inside the stress shadow. Then I perform a spatiotemporal search of PAS using the power-law fit methodology (i.e. c-value) and compare the optimal signal to the expected spatiotemporal extent of the theoretical signal. First I show that the optimal starting time and spatial extent of PAS are poorly controlled, due in part to the intrinsic properties of the c-value, but also to the random behavior of background seismicity. Second I show that theoretical PAS is identified by an optimal c-value (clear acceleration) only if the regional seismic activity (~ δ0b) is high and the noise ratio (δ−b/δ0b) is low. Otherwise the signal tends to disappear and the c-value becomes unstable. As a consequence, even if the power-law fit methodology is a simple approach to test the presence of PAS and can help provide a better understanding of the process engaged, it seems inadequate for robust systematic prospective forecasts. © 2008 Elsevier B.V. All rights reserved.
1. Introduction Accelerating Moment Release (AMR) has been proposed as a possible precursor to large earthquakes and is typically modeled by a simple power-law time-to-failure equation. Following Bufe and Varnes (1993), the relation is eðt Þ ¼ A þ Bðtf t Þm where tf is the time of the mainshock, B is negative and m is usually comprised between 0.3 and 0.7 (e.g. Mignan, 2007). For a convenient data analysis, ε(t) is chosen to be the cumulative Benioff strain but Mignan et al. (2007) show that the cumulative number of events should be preferred. The term Precursory Accelerating Seismicity (PAS) is used in the following study to distinguish the cumulative number of events k(t) from the cumulative Benioff strain ε(t) (AMR). It is important to note that the results obtained in this study for PAS
E-mail address: [email protected]. 0040-1951/$ – see front matter © 2008 Elsevier B.V. All rights reserved. doi:10.1016/j.tecto.2008.02.010
would be similar for AMR since the two power-laws have basically the same behavior. Despite the fact that AMR has been claimed to be observed in many regions in the world and for different tectonic settings (e.g. Brehm and Braile,1998; Jaume and Sykes,1999; Robinson, 2000; Mignan et al., 2006a; Papazachos et al., 2007), its existence is still debated. This can be explained by the fact that (1) AMR is not observed systematically before large earthquakes (e.g. Mignan et al. (2006a) failed to predict the 12 September 2007 M=8.4 Sumatra earthquake), (2) the duration and spatial extent of AMR are highly variable and do not seem to be linked to any characteristic of the subsequent mainshock, and finally (3) some spurious cases of AMR may be identified by the power-law fit methodology. The power-law fit methodology is used to quantify the degree of acceleration of precursory seismicity by comparing the root-meansquare error of a power-law fit to the root-mean-square error of a linear fit. This has been proposed by Bowman et al. (1998) and is defined as the curvature parameter C (or c-value) C ¼ powerlaw f it RMS error=linear f it RMS error
A. Mignan / Tectonophysics 452 (2008) 42–50
with the power exponent m required to be less than 1, C tends to zero if the seismicity is accelerating, otherwise C tends to 1 or more than 1 (i.e. seismicity linearly increasing or decelerating). The methodology is then applied for different spatial extents around the mainshock and for different time durations before the mainshock, in order to find the best AMR (or PAS). However, given sufficiently flexible search parameters, most optimization techniques will find some degree of acceleration in any data set. This bias can be quantified by using random synthetic seismicity catalogues (e.g. Bowman et al., 1998). Hardebeck et al. (submitted for publication) also examine whether or not AMR is a statistically significant precursor but conclude that apparent AMR results from a combination of data-fitting and the spatiotemporal clustering of earthquakes. However this only proves that clusters create artifacts (a cluster occurring at the end of the time series will lead to an apparent acceleration; a cluster occurring at the beginning of the time series will lead to an apparent deceleration). Mignan et al. (2006b) already discussed this issue and showed that declustering is necessary. This is in agreement with Mignan et al. (2007) who demonstrate that PAS (or AMR) corresponds in theory to a smooth acceleration through time and not to jumps of seismicity before the incoming mainshock. In the present article, I focus on the two other points which are at the origin of the skepticism around PAS existence (i.e. PAS not systematically observed, spatial extent and time duration highly variable). The originality of this work is that it studies theoretical synthetic catalogues where PAS exists, whereas other studies compare real seismicity catalogues, where PAS may or may not exist, with random synthetic catalogues, where there is no PAS. The purpose of this study is to determine, in the hypothesis that PAS exists, if it is likely to be fitted (using the c-value methodology). Theoretical catalogues are simulated based on some simple concepts first described by King and Bowman (2003), then developed by Mignan et al. (2007) and at present termed the Non-Critical Precursory Accelerating Seismicity Theory (NC PAST). 2. The Non-Critical Precursory Accelerating Seismicity Theory (NC PAST) 2.1. Previous formulations There is a general consensus in the scientific community to explain AMR/PAS by the theory of the Critical Point (e.g. Sammis and Sornette, 2002). Even studies that only search for AMR in seismicity catalogues assume that this precursor is due to critical processes (e.g. Bowman et al., 1998; Jaume and Sykes, 1999; Papazachos et al., 2007). However, another alternative has recently been proposed to explain the powerlaw behavior of AMR/PAS, based on the simple concept of elastic rebound (Reid, 1910). In that view, criticality is not necessary. The concept of elastic rebound in explaining AMR has first been proposed by King and Bowman (2003) in the Stress Accumulation model. Simulations show that AMR emerges from the background seismicity due to minor stress release as the whole region becomes sufficiently stressed for the mainshock to occur (i.e. stress loading on the fault, at depth). The Stress Accumulation model has been applied to real cases by Bowman and King (2001) and Mignan et al. (2006a,b). Mignan et al. (2007) developed a new approach, slightly different from the original Stress Accumulation model, to clearly identify the origin of the AMR power-law. Indeed, the origin of AMR in the Stress Accumulation model remains unclear because (1) simulations use complicated stress patterns (i.e. loading stress + tectonic stress memory, see details in King and Bowman (2003)) and (2) the role of the frequency–magnitude relation on the AMR power-law is unknown. As a consequence, the Stress Accumulation model is generally considered as a kind of a Critical Point process (Tiampo and Anghel, 2006), equivalent to a percolation model (Sammis and
43
Sornette, 2002). Mignan et al. (2007) demonstrated that criticality is not necessary and that the power-law behavior of PAS (and consequently of AMR) is purely geometrical. This article is the continuation of this new approach, which I term Non-Critical Precursory Accelerating Seismicity Theory (or NC PAST). 2.2. New formulation In the NC PAST, PAS can emerge from the background seismicity because of the decrease, due to loading, of the size of a stress shadow due to a previous earthquake. In other words, PAS is due to the increase of the size of the region where background seismicity is allowed to occur. Fig. 1 illustrates the spatiotemporal distribution of background seismicity in the stress field of a loading fault, represented by a source point at r = 0. At t = 0, an earthquake occurs on the fault, which creates a stress shadow. The size of the stress shadow decreases through time to finally disappear at t = tf allowing a new earthquake to occur on the fault (i.e. mainshock). The spatiotemporal evolution of the stress field is controlled by a constant stress loading on the fault during its seismic cycle. For display reasons, the spatial evolution of the stress shadow is represented for only one dimension of space, r (see Fig. 3 for a two-dimensional representation in space of the stress shadow). Following Mignan et al. (2007), the contour limit of the stress shadow (i.e. the size of the stress shadow) is defined by the equation r⁎(t, σ = σ⁎) where σ⁎ approaches the failure stress σf and means that events are more likely to occur where r ≥ r⁎. The relation is of the form r⁎(t) ∞ (tf − t)1/3. The density per unit of space and time of events that compose the background seismicity is δ0b outside the stress shadow and δ−b inside the stress shadow, with δ0b N N δ−b. Mignan et al. (2007) demonstrated analytically that for a fixed region in space (optimum from r = 0 to r~r⁎max), the cumulative number of events k(t) that compose the background seismicity (δ0b) increases as a power-law function though time before the mainshock, as k(t) ∞ −(tf − t)D/3 with D being the fractal dimension of the associated fault network. Events that occur in the stress shadow (δ−b) tend to hide the precursory seismicity pattern (i.e. background seismicity noise). Fig. 2 represents the relation between the density of background seismicity δb and the stress σ. In an ideal case (left panel), events would only occur where the stress is positive (i.e. aftershocks) and regions of negative stress would remain inactive (i.e. quiescence). In a more realistic view (right panel), events are not influenced by very small stress changes (near-zero values) and occur randomly in space
Fig. 1. Spatiotemporal distribution of background seismicity in the stress field of a loading fault (represented by a source point, at r = 0). At t = 0, an earthquake occurs on the fault, which creates a stress shadow. The size of the stress shadow decreases thought time to finally disappear at t = tf allowing a new earthquake to occur on the fault. For display reasons, the spatial evolution of the stress shadow is represented for only one dimension of space, r. Following Mignan et al. (2007), the contour limit of the stress shadow (i.e. the size of the stress shadow) is defined by the equation r⁎(t, σ = σ⁎) where σ⁎ approaches the failure stress σf and means that events are more likely to occur where r ≥ r⁎. The density per unit of space and time of events that compose the background seismicity is δ0b outside the stress shadow and δ−b inside the stress shadow, with δ0b N N δ−b.
44
A. Mignan / Tectonophysics 452 (2008) 42–50
Fig. 2. Schematic representation of the relation between the density of background seismicity δb and the stress σ. In an ideal case (left), events would only occur where the stress is positive (i.e. aftershocks) and regions of negative stress would remain inactive (i.e. quiescence). In a more realistic view (right), events are not influenced by very small stress changes (near-zero values) and occur randomly in space and time. These events compose the background seismicity (of density δ0b), which is generally observed in far field. Aftershocks (of density δ+b) occur where the stress is sufficient to trigger new events. In lobes of negative stress (i.e. stress shadow), there is a relative quiescence (of density δ−b) with only a few events being able to occur, such as δ+b N Nδ0b N N δ−b. The study of PAS only involves the spatiotemporal evolution of the stress limit σ⁎, which separates the background seismicity from erratic events in the stress shadow.
and time. These events compose the background seismicity of density δ0b, which is generally observed in far field. Aftershocks occur where the stress is sufficient to trigger new events. This triggering process increases the density of events to δ+b in lobes of positive stress. In lobes of negative stress (i.e. stress shadow), there is a relative quiescence with only a few events being able to occur for several reasons (e.g. fast stressing rate, transient effects, complex stress history). However, the density δ−b is supposed to be low with δ+b N N δ0b N N δ−b. Note that there is no reason for the density of events δ to evolve by step at given stress values. In particular for aftershocks, δ should increase continuously with increasing stress (I take the simple case where only static stress changes occur). However, this assumption has no consequences here since the study of PAS only involves the spatiotemporal evolution of the stress limit σ⁎, which separates the background seismicity (random in space) from erratic events in the stress shadow (seismicity noise that is also assumed random in space). I assume in this study that PAS (if it exists) is explained by the NC PAST. Other possible explanations (e.g. based on the theory of the Critical Point) are dismissed for the following reasons: (1) The NC PAST seems the most simple way to explain the powerlaw time-to-failure equation. It is based on simple and recognized concepts, such as elastic rebound and static stress transfer, and the origin of the power-law is purely geometrical. Therefore, I follow Ockham's Razor from which, of two equivalent theories or explanations, all other things being equal, the simpler one is to be preferred. (2) The NC PAST explains, in addition to PAS, other proposed seismicity precursors such as quiescence and the Mogi “doughnut” (e.g. Mogi, 1969). Note that this is already proposed in the Stress Accumulation model (King and Bowman, 2003). (3) Contrary to phenomena due to critical processes, which are difficult to catch in a complex dynamic system, PAS can be simply captured in space and time in the NC PAST view. This permits to determine how this precursor can be extracted, and thus to improve forecasting techniques. This aspect is used in the following simulations to test the power-law fit methodology (i.e. c-value). 3. Theoretical synthetic seismicity catalogues 3.1. Simulation methodology Theoretical synthetic catalogues are created based on the principles of the NC PAST. The background seismicity density δ0b and the noise ratio δ−b/δ0b are adjustable parameters that allow to test the consistency of PAS, depending respectively on the regional seismic
activity and on the effect of the stress shadow on seismicity. Coulomb stress calculations are realized using the source code AlmondX (courtesy of G. C. P. King), based on the dislocation theory in an elastic half-space (Okada, 1992). Events location x,y and time t are controlled by the stress contour limit σ⁎(x,y,t). The stress field σ(x,y,t) is calculated for a dip-slip fault (black segment) located in the center of a grid [0,1; 0,1]. The time interval is [t0 = 0; tf = 20]. For each time step ti (fixed to 1 year), events are random uniform in space and time (in this year interval) with a density δ0b outside the stress contour σ⁎(x,y,ti) and δ−b inside the stress contour σ⁎(x,y,ti) (Fig. 3, left panel). A more realistic representation of the behavior of background seismicity would take into account spatiotemporal clustering (i.e. ETAS model (e.g. Ogata and Zhuang, 2006) instead of a simple random uniform event distribution). However, since seismicity must be declustered in order to apply the power-law fit methodology, the latter approach seems more appropriate and straightforward. The stress field σ(x,y,ti) results from a displacement on the fault (imax − i)ΔD with ΔD being a constant displacement increment. At t0, an earthquake occurs on the fault with a displacement imaxΔD producing a stress shadow, which then decreases through time because of the component − iΔD (Fig. 3, right). This negative component corresponds to the linear superposition of a loading stress field by following the concept of back-slip model (Savage, 1983). The value of σ⁎ is fixed such that the last contour limit σ⁎(x,y,tf) is close to the fault segment (bL/10 with L the fault length), meaning that the fault is ready to fail (i.e. earthquake on the fault at t ~ tf, represented by a star on Fig. 3). Note that the value of σ⁎ depends on the fault mechanism, the displacement increment ΔD and the number of time steps ti. The relation r⁎(t, σ = σ⁎) that constrains the size of the stress shadow through time in one dimension (Mignan et al., 2007) can be applied to a two-dimensional case with r⁎(t, σ = σ⁎) valid along the gradient vectors of the stress field (two gradients represented on Fig. 3, right panel). However simulations are preferred since no simple analytical formulation can be given here. Fig. 4 shows an example of accelerating seismicity pattern (C = 0.29) observed in a theoretical synthetic catalogue. Events that compose the precursory signal are located in the region inside the stress contour σ⁎(t = t0), which corresponds to the maximum extent of the stress shadow, produced by the fault rupture at t = t0. The cumulative number of events λ increases as a power-law function through time, in agreement with Mignan et al. (2007). 3.2. Power-law fit methodology To detect accelerating seismicity patterns in a data set (here in theoretical synthetic catalogues), it is necessary to isolate the precursory events in a space/time window (also in a magnitude window in real cases, but not required here). Many studies use simple shapes such as circles (e.g. Brehm and Braile, 1998; Bowman et al., 1998; Robinson, 2000; Zoller et al., 2001) or ellipses (e.g. Papazachos et al., 2007) and a few other studies, based on the Stress Accumulation model (~NC PAST), use stress contours to define the region of interest (e.g. Mignan et al., 2006b; Bowman and King, 2001). In theoretical synthetic catalogues, PAS is located, by definition of the NC PAST, in the region delimited by a stress shadow. However in this study, since it is clearly impossible to determine the exact region of interest in real cases (i.e. stress contour depending on fault geometry and mechanism, on interactions with other active structures, etc…), PAS is identified in circular regions. This allows to obtain more realistic results from the optimization procedure. Fig. 5 (left panel) illustrates the optimization procedure employed for the systematic search of accelerating seismicity patterns in space and time. Patterns are searched for in concentric circles of different radii R centered on the fault and for different starting times tstart. The final time is fixed to tf. For each step (10 in space with 0.025 b R b 0.250
A. Mignan / Tectonophysics 452 (2008) 42–50
45
Fig. 3. Simulation of theoretical synthetic seismicity catalogues (NC PAST based). Events location x,y and time t are controlled by the stress contour limit σ⁎(x,y,t). The stress field σ(x,y,t) is calculated for a dip-slip fault (black segment) located in the center of a grid [0,1; 0,1]. The time interval is [t0 =0; tf = 20]. For each time step ti (fixed to 1 year), events are random uniform in space and time (in this year interval) with a density δ0b outside the stress contour σ⁎(x,y,ti) and δ−b inside the stress contour σ⁎(x,y,ti) (left panel). The value of σ⁎ is fixed such as the last contour limit σ⁎(x,y,tf) is close to the fault segment (bL/10 with L the fault length), meaning that the fault is ready to fail (i.e. earthquake on the fault at t ~tf, represented by a star). The relation r⁎(t, σ =σ⁎) that constrains the size of the stress shadow through time in one dimension (Mignan et al., 2007) can be applied to the two-dimensional case with r⁎(t, σ =σ⁎) being valid along the gradient vectors of the stress field (two gradients represented, right panel).
and 10 in time with 0 b tstart b 18), the c-value is calculated and the result plotted in the spatiotemporal search matrix or c-value matrix (right panel). In previous studies using the same optimization procedure (e.g. Sammis et al., 2004; Mignan et al., 2006a), the accelerating seismicity pattern is chosen as the signal with the lowest c-value (note that signals with a low c-value but composed of less than 10–15 events are generally dismissed). In theoretical synthetic catalogues, the optimal spatial extent and time duration of PAS are already known. The optimal starting time is t0 = 0 and in the example of Fig. 5, PAS is located in a circle of radius 0.10 b R b 0.15 (bold circles). R can only be approximated since no circle can fit the exact location of PAS (dark grey region). 4. Results
time in the interval 0 b t0 b10 (in theory, 0 b t0 b 3) and in a circle with 0.075 b R b 0.175 (in theory, 0.10 b R b 0.15). In addition, it must be noted that some instabilities can appear with a low c-value observed in more than one region (i.e. several blue patches). In such a case, it becomes difficult to determine the correct location of PAS. Note that Figs. A1–A5 clearly show the extreme variability of the c-value for different simulations with the same input parameters δ0b and δ−b/δ0b. The spatiotemporal variability/instability of the c-value can be explained by the fact that the power-law fit is very sensitive to the few events that occur in the first part of the acceleration (see Fig. 4). Because of the random behavior of events in space and time in the background, the starting point of the power-law curve is poorly controlled, leading to the variations observed among the c-value matrices. This tends to prove that the c-value does not permit to clearly determine the duration and spatial extent of PAS. This is due in part to the
4.1. Weakness of the power-law fit methodology Fig. 5 (right panel) shows an example of c-value matrix, which corresponds to a stack of 50 simulation results. In that particular case, the best acceleration (C b 0.3, blue region) is found for a starting time 0 b t0 b 3 and in a circle with a radius R comprised between 0.10 and 0.15 (bold circles on the left panel), which corresponds perfectly to the spatiotemporal window where the theoretical acceleration is located. This result is obtained thanks to the stacking method, which allows to increase the signal-to-noise ratio of the analysis, as proposed in some other forecasting techniques, such as Pattern Informatics (Holliday et al., 2005; Tiampo et al., 2002). Here the stacking method consists of creating a c-value matrix in which each value C(R, t0) corresponds to the mean of the values C(R, t0) computed from n simulations. The c-value matrices obtained for the 50 simulations (δ0b = 1000, δ−b/ δ0b = 0, see next section for different input parameters) are represented on Fig. A1 and four examples are represented on Fig. 6. Although PAS remains the same in theory (same input parameters δ0b and δ−b/δ0b), the random behavior of background seismicity that composes the precursory signal makes each pattern unique (example of Fig. 4). Differences do not change the shape of the acceleration at the first order but are sufficient to make significant changes on the power-fit result. On Fig. 6, the lowest c-value (C b 0.3) is obtained for a starting
Fig. 4. Example of accelerating seismicity pattern (C = 0.29) observed in a theoretical synthetic catalogue. Events that compose the precursory signal are located in the region inside the stress contour σ⁎(t = t0), which corresponds to the maximum extent of the stress shadow, produced by the fault rupture at t = t0. The cumulative number of events λ increases as a power-law function through time, in agreement with Mignan et al. (2007). The power-law fit and linear fit used to calculate the c-value are represented.
46
A. Mignan / Tectonophysics 452 (2008) 42–50
Fig. 5. Method employed for a systematic search of accelerating seismicity patterns in space and time. Left panel: Patterns are searched in circles of different radii centered on the fault and for different starting times. The final time is fixed to tf. Right panel: Spatiotemporal search matrix or c-value matrix. C-values are computed for 10 starting times in the interval [0; 18] and for 10 radii in the interval [0.025; 0.250]. In this particular example, the result corresponds to a stack of c-value matrices obtained from 50 simulations (δ0b = 1000, δ−b/δ0b =0.).
inherent weakness of the power-law fit but also to the random behavior of events that compose PAS. If the first issue (linked to the fit methodology) could only be resolved by adopting a new approach to extract PAS (not the purpose of this article), the second issue (linked to the random behavior of events) can be resolved by stacking many data sets. 4.2. Consistency of PAS As already discussed before, the existence of AMR (~PAS) is usually questioned on the basis that the duration and spatial extent of AMR (~PAS)
are highly variable and do not seem to be linked to any characteristic of the subsequent mainshock. This argument seems at present irrelevant since this variability is also observed in theoretical synthetic catalogues where PAS is known to exist. The other argument against the existence of AMR (~PAS) is that it is not observed systematically before large earthquakes. This argument cannot be refuted by the results of the previous simulations (δ0b =1000, δ−b/δ0b =0.), where a clear acceleration (i.e. low c-value) is found in the majority of cases (see Figs. 6 and A1). Therefore some new simulations are realized with different input parameters (δ0b and δ−b/δ0b) in order to test the consistency of PAS. First,
Fig. 6. Four examples of c-value matrix. Although the input parameters of the different simulations are the same (i.e. δ0b = 1000, δ−b/δ0b = 0), the power-law fit result is highly variable. This is due in part to the inherent weakness of the power-law fit methodology but also to the random behavior of background seismicity that composes PAS. See Fig. A1 for the c-value matrices of 50 simulations.
A. Mignan / Tectonophysics 452 (2008) 42–50 Fig. 7. Role of the background seismicity density δ0b in the formation of PAS (noise ratio δ−b/δ0b fixed to zero). From left to right, δ0b = 1000, δ0b = 500 and δ0b = 100 (i.e. number of events in the unit square [0,1; 0,1], per year). For each value of δ0b, the upper panel shows an example of theoretical PAS. The central panel shows the probability density function (PDF) of the corresponding c-value, for 100 simulations. The lower panel shows the c-value matrix obtained after stacking the results of 50 simulations (All c-value matrices are represented in Figs. A1–A3 respectively). 47
48 A. Mignan / Tectonophysics 452 (2008) 42–50 Fig. 8. Role of the noise ratio δ−b/δ0b in the formation of PAS (background density δ0b fixed to 1000). From left to right, δ−b/δ0b = 0.0, δ−b/δ0b = 0.2 and δ−b/δ0b = 0.5. For each value of δ−b/δ0b, the upper panel shows an example of theoretical PAS. The central panel shows the probability density function (PDF) of the corresponding c-value, for 100 simulations. The lower panel shows the c-value matrix obtained after stacking the results of 50 simulations (All c-value matrices are represented in Figs. A1, A4 and A5 respectively).
A. Mignan / Tectonophysics 452 (2008) 42–50
the possible role of a decrease of the regional seismic activity on PAS is examined by changing the value of δ0b (noise ratio δ−b/δ0b fixed to zero). Fig. 7 illustrates the evolution of PAS for δ0b = 1000, 500 and 100 (i.e. number of events in the unit square [0,1; 0,1], per year). For each value of δ0b, the upper panel shows an example of theoretical PAS, which is defined from its exact spatiotemporal extent. The number of events in the signal is λf. The central panel shows the probability density function (PDF) of the corresponding c-value, for 100 simulations. The lower panel shows the c-value matrix obtained after stacking the results of 50 simulations (all c-value matrices are represented in Figs. A1–A3 respectively). It can be observed that a change of the background seismicity density δ0b has a clear impact on the consistency of PAS. The lower the density δ0b, the more the precursory signal tends to disappear. On the c-value PDF (central panel), the coefficient of variation increases and the mean value shifts to higher c-values. This phenomenon is not a consequence of the weakness of the power-law fit since it can still be observed after stacking many c-value matrices (lower panel) but is due to the way PAS forms in the NC PAST. In that view, an ideal acceleration would form if δ0b tends to infinity. Mignan et al. (2007) defined the power-law time-to-failure equation for a continuous quantity, linked to the continuous decrease of the size of a stress shadow. Since the precursory signal is composed of punctual events, some information is irremediably lost. With δ0b decreasing, the signal fluctuates more and more around the ideal power-law curve, as observed on the upper panel of Fig. 7. Fig. 8 shows the impact of an increase of the noise ratio δ−b/δ0b on PAS for δ−b/δ0b = 0.0, 0.2 and 0.5 (with δ0b fixed to 1000). In the same way, the higher the noise ratio δ−b/δ0b, the more the precursory signal tends to disappear. In that case, the theoretical PAS is still present (clear for δ0b = 1000) but a linear trend is superimposed to the power-law. Therefore the power-law fit methodology is less efficient. With δ−b/ δ0b = 0.2, the best c-value obtained in the spatiotemporal search matrix (lower panel) increases from 0.3 to 0.5 and with δ−b/δ0b = 0.5, the signal becomes very weak (C ~ 0.7), which would not be considered as a reliable signal in a real case. These two tests permit to explain why PAS (or AMR) is not systematically observed before large earthquakes. If the NC PAST is correct, PAS would only be clearly observed in regions of high seismic activity. Moreover, quiescence phenomena are not systematically observed in stress shadows. Hardebeck et al. (1998) showed that only 65% of aftershocks occurred in regions of stress increase related to the 1992 Landers earthquake, 60% to the 1994 Northridge earthquake. Therefore between 35 and 40% of all aftershocks occurred in regions of decreased stress, part of these aftershocks certainly occurred in the stress shadow delimited by the negative stress value σ⁎. This tends to prove that the noise ratio δ−b/δ0b might be high in real seismicity catalogues, leading to precursory patterns difficult to catch. 5. Discussion The purpose of this article has not been to prove or to verify that PAS exists but to show, assuming that the NC PAST is correct, that one can explain why PAS is not observed systematically before large earthquakes and why the spatial extent and duration of the signal seem unstable using the c-value. It is important to note that Bowman et al. (1998) clearly formalized the AMR concept into a testable hypothesis thanks to the c-value (noted by Hardebeck et al., submitted for publication). The c-value provides a simple, computationally inexpensive parameter to search for AMR/PAS patterns and thus has been used in many studies since Bowman et al. (1998) (see review in Mignan, 2007). Nowadays, authors rarely make the distinction between c-value and AMR, leading irremediably to some misconceptions. I have shown in this study that an unstable c-value does not mean that AMR/PAS does not exist. I have demonstrated why the signal can
49
be weak and what parameters control its appearance (i.e. the background seismic activity δ0b and the background seismicity noise ratio δ−b/δ0b). These results lead to three important conclusions: (1) The c-value is inadequate for robust systematic prospective forecasts. The c-value only gives consistent results for regions of high background seismic activity (δ0b high) where stress shadows are clearly defined (δ−b/δ0b low). This case is rarely encountered in real seismicity catalogues. It is interesting to note that stacking c-value matrices permits to improve the power-law fit results. However this can only be done in retrospect and requires that PAS patterns due to different large earthquakes have the same spatiotemporal extent, which is also unlikely. (2) PAS/AMR possibly exists if the NC PAST is correct. This new theory makes specific predictions for the spatiotemporal extent of PAS and explains why the signal can appear or not. Future studies based on the NC PAST should focus on some new methodologies that would permit to extract PAS in a better way. The c-value has been proved to be simple and useful in the past but new results show that more complicated algorithms are necessary to possibly forecast large earthquakes. (3) Another alternative would be to improve the seismic network detection threshold. The NC PAST shows indeed that it is the number of events which is important and not the magnitude for the appearance of PAS. This could be done by increasing the density of simple one component narrow band seismometers in shallow boreholes. At present the trend in networks is to have fewer 3 component instruments with the result that detection threshold has declined in many places such as South California and Greece (G. C. P. King, pers. comm.). Acknowledgements The author acknowledges helpful reviews by Kristy Tiampo and Geoffrey King as well as Editor Tom Parsons. The author also thanks Geoffrey King for the source code AlmondX used for the Coulomb stress calculations. Appendix A. Supplementary data Supplementary data associated with this article can be found, in the online version, at doi:10.1016/j.tecto.2008.02.010. References Bowman, D.D., King, G.C.P., 2001. Accelerating seismicity and stress accumulation before large earthquakes. Geophys. Res. Lett. 28, 4039–4042. doi:10.1029/2001GL013022.2001. Bowman, D.D., Ouillon, G., Sammis, C.G., Sornette, A., Sornette, D., 1998. An observational test of the critical earthquake concept. J. Geophys. Res. 103, 24,359–24,372. Brehm, D.J., Braile, L.W., 1998. Intermediate-term earthquake prediction using precursory events in the New Madrid Seismic Zone. Bull. Seismol. Soc. Am. 88, 564–580. Bufe, C.G., Varnes, D.J., 1993. Predictive modeling of the seismic cycle of the greater San Francisco Bat region. J. Geophys. Res. 98, 9871–9883. Hardebeck, J.L., Nazareth, J.J., Hauksson, E., 1998. The static stress change triggering model: constraints from two southern California aftershock sequences. J. Geophys. Res. 103, 24,427–24,437. Hardebeck, J.L., Felzer, K.R., Michael, A.J. (submitted for publication), Rigorous observational tests contradict the accelerating moment release hypothesis. Holliday, J.R., Nanjo, K.Z., Tiampo, K.F., Rundle, J.B., Turcotte, D.L., 2005. Earthquake forecasting and its verification. Nonlinear Process Geophys. 12, 965–977. Jaume, S.C., Sykes, L.R.,1999. Evolving towards a critical point: a review of accelerating seismic moment/energy release prior to large and great earthquakes. Pure Appl. Geophys. 155, 279–306. King, G.C.P., Bowman, D.D., 2003. The evolution of regional seismicity between large earthquakes. J. Geophys. Res. 108. doi:10.1029/2001JB000783. Mignan, A., 2007. The Stress Accumulation model: accelerating moment release & seismic hazard. Adv. Geophys. 49. doi:10.1016/S0065-2687(07)49002-1. Mignan, A., King, G.C.P., Bowman, D.D., Lacassin, R., Dmowska, R., 2006a. Seismic activity in the Sumatra-Java region prior to the December 26, 2004 (Mw = 9.0–9.3) and March 28, 2005 (Mw = 8.7) earthquakes. Earth Planet. Sci. Lett. 244, 639–654.
50
A. Mignan / Tectonophysics 452 (2008) 42–50
Mignan, A., Bowman, D.D., King, G.C.P., 2006b. An observational test of the origin of accelerating moment release before large earthquakes. J. Geophys. Res. 111. doi:10.1029/2006JB004374. Mignan, A., King, G.C.P., Bowman, D.D., 2007. A mathematical formulation of accelerating moment release based on the stress accumulation model. J. Geophys. Res. 112. doi:10.1029/2006JB004671. Mogi, K., 1969. Some features of recent seismic activity in and near Japan (2), Activity before and after large earthquakes. Bull. Earthq. Res. Inst. Univ. Tokyo 47, 395–417. Ogata, Y., Zhuang, J., 2006. Space–time ETAS models and an improved extension. Tectonophysics 413, 13–23. Okada, Y., 1992. Internal deformation due to shear and tensile faults in a half-space. Bull. Seismol. Soc. Am. 82, 1018–1040. Papazachos, B.C., Karakaisis, G.F., Papazachos, C.B., Scordilis, E.M., 2007. Evaluation of the results for an intermediate-term prediction of the 8 January 2006 Mw = 6.9 Cythera earthquake in the Southwestern Aegean. Bull. Seismol. Soc. Am. 97, 347–352. doi:10.1785/0120060075. Reid, H.F., 1910. The mechanics of the earthquake, the California earthquake of April 18, 1906. Report State Investig. Comm., vol. 2. Carnegie Inst., Washington.
Robinson, R., 2000. A test of the precursory Accelerating Moment Release model on some recent New Zealand earthquakes. Geophys. J. Int. 140, 568–576. Sammis, C.G., Sornette, D., 2002. Positive feedback, memory, and the predictability of earthquakes. Colloq. Proc. - Natl. Acad. Sci. 99, 2501–2508. Sammis, C.G., Bowman, D.D., King, G.C.P., 2004. Anomalous seismicity and accelerating moment release preceding the 2001 and 2002 earthquakes in northern Baja California, Mexico. Pure Appl. Geophys. 161, 2369–2378. Savage, J.C., 1983. A dislocation model of strain accumulation and release at a subduction zone. J. Geophys. Res. 88, 4984–4996. Tiampo, K.F., Anghel, M., 2006. Introduction to special issue: critical point theory and space–time pattern formation in precursory seismicity. Tectonophysics 413, 1–3. Tiampo, K.F., Rundle, J.B., McGinnis, S., Klein, W., 2002. Pattern dynamics and forecast methods in seismically active regions. Pure Appl. Geophys. 159, 2429–2467. Zoller, G., Hainzl, S., Kurths, J., 2001. Observation of growing correlation length as an indicator for critical point behavior prior to large earthquakes. J. Geophys. Res. 106, 2167–2176.
Contents lists available at ScienceDirect
Tectonophysics j o u r n a l h o m e p a g e : w w w. e l s ev i e r. c o m / l o c a t e / t e c t o
Non-Critical Precursory Accelerating Seismicity Theory (NC PAST) and limits of the power-law fit methodology A. Mignan Risk Management Solutions, Science and Technology Research, Peninsular House, 30 Monument Street, London EC3R 8NB, United Kingdom
a r t i c l e
i n f o
Article history: Received 18 September 2007 Received in revised form 2 February 2008 Accepted 13 February 2008 Available online 4 March 2008 Keywords: Earthquake Forecast Precursory seismicity AMR Stress loading Stress shadow Synthetic catalogue
a b s t r a c t The hypothesis that Accelerating Moment Release (AMR) is a precursor to large earthquakes is still debated. On one hand, AMR has been claimed to be observed in many cases and on the other hand, it has been proposed that apparent AMR is only due to data-fitting. The debate is in general focused on the validity of the c-value (curvature parameter), which permits to quantify AMR (i.e. cumulative Benioff strain through time), or more generally Precursory Accelerating Seismicity (PAS, i.e. cumulative number of events through time). Contrary to previous studies, which compare c-value optimization in real seismicity catalogues and in random synthetic catalogues, I test c-value optimization in theoretical synthetic catalogues. In that particular case, I assume that PAS exists and that it can be explained by the Non-Critical Precursory Accelerating Seismicity Theory (NC PAST). This theory demonstrates that PAS can emerge from the background seismicity because of the decrease, due to loading, of the size of a stress shadow due to a previous earthquake. I improve the NC PAST by integrating the following characteristics of the background seismicity, (1) the density of random events outside the stress shadow δ0b and (2) the noise ratio δ−b/δ0b, with δ−b being the density of random events inside the stress shadow. Then I perform a spatiotemporal search of PAS using the power-law fit methodology (i.e. c-value) and compare the optimal signal to the expected spatiotemporal extent of the theoretical signal. First I show that the optimal starting time and spatial extent of PAS are poorly controlled, due in part to the intrinsic properties of the c-value, but also to the random behavior of background seismicity. Second I show that theoretical PAS is identified by an optimal c-value (clear acceleration) only if the regional seismic activity (~ δ0b) is high and the noise ratio (δ−b/δ0b) is low. Otherwise the signal tends to disappear and the c-value becomes unstable. As a consequence, even if the power-law fit methodology is a simple approach to test the presence of PAS and can help provide a better understanding of the process engaged, it seems inadequate for robust systematic prospective forecasts. © 2008 Elsevier B.V. All rights reserved.
1. Introduction Accelerating Moment Release (AMR) has been proposed as a possible precursor to large earthquakes and is typically modeled by a simple power-law time-to-failure equation. Following Bufe and Varnes (1993), the relation is eðt Þ ¼ A þ Bðtf t Þm where tf is the time of the mainshock, B is negative and m is usually comprised between 0.3 and 0.7 (e.g. Mignan, 2007). For a convenient data analysis, ε(t) is chosen to be the cumulative Benioff strain but Mignan et al. (2007) show that the cumulative number of events should be preferred. The term Precursory Accelerating Seismicity (PAS) is used in the following study to distinguish the cumulative number of events k(t) from the cumulative Benioff strain ε(t) (AMR). It is important to note that the results obtained in this study for PAS
E-mail address: [email protected]. 0040-1951/$ – see front matter © 2008 Elsevier B.V. All rights reserved. doi:10.1016/j.tecto.2008.02.010
would be similar for AMR since the two power-laws have basically the same behavior. Despite the fact that AMR has been claimed to be observed in many regions in the world and for different tectonic settings (e.g. Brehm and Braile,1998; Jaume and Sykes,1999; Robinson, 2000; Mignan et al., 2006a; Papazachos et al., 2007), its existence is still debated. This can be explained by the fact that (1) AMR is not observed systematically before large earthquakes (e.g. Mignan et al. (2006a) failed to predict the 12 September 2007 M=8.4 Sumatra earthquake), (2) the duration and spatial extent of AMR are highly variable and do not seem to be linked to any characteristic of the subsequent mainshock, and finally (3) some spurious cases of AMR may be identified by the power-law fit methodology. The power-law fit methodology is used to quantify the degree of acceleration of precursory seismicity by comparing the root-meansquare error of a power-law fit to the root-mean-square error of a linear fit. This has been proposed by Bowman et al. (1998) and is defined as the curvature parameter C (or c-value) C ¼ powerlaw f it RMS error=linear f it RMS error
A. Mignan / Tectonophysics 452 (2008) 42–50
with the power exponent m required to be less than 1, C tends to zero if the seismicity is accelerating, otherwise C tends to 1 or more than 1 (i.e. seismicity linearly increasing or decelerating). The methodology is then applied for different spatial extents around the mainshock and for different time durations before the mainshock, in order to find the best AMR (or PAS). However, given sufficiently flexible search parameters, most optimization techniques will find some degree of acceleration in any data set. This bias can be quantified by using random synthetic seismicity catalogues (e.g. Bowman et al., 1998). Hardebeck et al. (submitted for publication) also examine whether or not AMR is a statistically significant precursor but conclude that apparent AMR results from a combination of data-fitting and the spatiotemporal clustering of earthquakes. However this only proves that clusters create artifacts (a cluster occurring at the end of the time series will lead to an apparent acceleration; a cluster occurring at the beginning of the time series will lead to an apparent deceleration). Mignan et al. (2006b) already discussed this issue and showed that declustering is necessary. This is in agreement with Mignan et al. (2007) who demonstrate that PAS (or AMR) corresponds in theory to a smooth acceleration through time and not to jumps of seismicity before the incoming mainshock. In the present article, I focus on the two other points which are at the origin of the skepticism around PAS existence (i.e. PAS not systematically observed, spatial extent and time duration highly variable). The originality of this work is that it studies theoretical synthetic catalogues where PAS exists, whereas other studies compare real seismicity catalogues, where PAS may or may not exist, with random synthetic catalogues, where there is no PAS. The purpose of this study is to determine, in the hypothesis that PAS exists, if it is likely to be fitted (using the c-value methodology). Theoretical catalogues are simulated based on some simple concepts first described by King and Bowman (2003), then developed by Mignan et al. (2007) and at present termed the Non-Critical Precursory Accelerating Seismicity Theory (NC PAST). 2. The Non-Critical Precursory Accelerating Seismicity Theory (NC PAST) 2.1. Previous formulations There is a general consensus in the scientific community to explain AMR/PAS by the theory of the Critical Point (e.g. Sammis and Sornette, 2002). Even studies that only search for AMR in seismicity catalogues assume that this precursor is due to critical processes (e.g. Bowman et al., 1998; Jaume and Sykes, 1999; Papazachos et al., 2007). However, another alternative has recently been proposed to explain the powerlaw behavior of AMR/PAS, based on the simple concept of elastic rebound (Reid, 1910). In that view, criticality is not necessary. The concept of elastic rebound in explaining AMR has first been proposed by King and Bowman (2003) in the Stress Accumulation model. Simulations show that AMR emerges from the background seismicity due to minor stress release as the whole region becomes sufficiently stressed for the mainshock to occur (i.e. stress loading on the fault, at depth). The Stress Accumulation model has been applied to real cases by Bowman and King (2001) and Mignan et al. (2006a,b). Mignan et al. (2007) developed a new approach, slightly different from the original Stress Accumulation model, to clearly identify the origin of the AMR power-law. Indeed, the origin of AMR in the Stress Accumulation model remains unclear because (1) simulations use complicated stress patterns (i.e. loading stress + tectonic stress memory, see details in King and Bowman (2003)) and (2) the role of the frequency–magnitude relation on the AMR power-law is unknown. As a consequence, the Stress Accumulation model is generally considered as a kind of a Critical Point process (Tiampo and Anghel, 2006), equivalent to a percolation model (Sammis and
43
Sornette, 2002). Mignan et al. (2007) demonstrated that criticality is not necessary and that the power-law behavior of PAS (and consequently of AMR) is purely geometrical. This article is the continuation of this new approach, which I term Non-Critical Precursory Accelerating Seismicity Theory (or NC PAST). 2.2. New formulation In the NC PAST, PAS can emerge from the background seismicity because of the decrease, due to loading, of the size of a stress shadow due to a previous earthquake. In other words, PAS is due to the increase of the size of the region where background seismicity is allowed to occur. Fig. 1 illustrates the spatiotemporal distribution of background seismicity in the stress field of a loading fault, represented by a source point at r = 0. At t = 0, an earthquake occurs on the fault, which creates a stress shadow. The size of the stress shadow decreases through time to finally disappear at t = tf allowing a new earthquake to occur on the fault (i.e. mainshock). The spatiotemporal evolution of the stress field is controlled by a constant stress loading on the fault during its seismic cycle. For display reasons, the spatial evolution of the stress shadow is represented for only one dimension of space, r (see Fig. 3 for a two-dimensional representation in space of the stress shadow). Following Mignan et al. (2007), the contour limit of the stress shadow (i.e. the size of the stress shadow) is defined by the equation r⁎(t, σ = σ⁎) where σ⁎ approaches the failure stress σf and means that events are more likely to occur where r ≥ r⁎. The relation is of the form r⁎(t) ∞ (tf − t)1/3. The density per unit of space and time of events that compose the background seismicity is δ0b outside the stress shadow and δ−b inside the stress shadow, with δ0b N N δ−b. Mignan et al. (2007) demonstrated analytically that for a fixed region in space (optimum from r = 0 to r~r⁎max), the cumulative number of events k(t) that compose the background seismicity (δ0b) increases as a power-law function though time before the mainshock, as k(t) ∞ −(tf − t)D/3 with D being the fractal dimension of the associated fault network. Events that occur in the stress shadow (δ−b) tend to hide the precursory seismicity pattern (i.e. background seismicity noise). Fig. 2 represents the relation between the density of background seismicity δb and the stress σ. In an ideal case (left panel), events would only occur where the stress is positive (i.e. aftershocks) and regions of negative stress would remain inactive (i.e. quiescence). In a more realistic view (right panel), events are not influenced by very small stress changes (near-zero values) and occur randomly in space
Fig. 1. Spatiotemporal distribution of background seismicity in the stress field of a loading fault (represented by a source point, at r = 0). At t = 0, an earthquake occurs on the fault, which creates a stress shadow. The size of the stress shadow decreases thought time to finally disappear at t = tf allowing a new earthquake to occur on the fault. For display reasons, the spatial evolution of the stress shadow is represented for only one dimension of space, r. Following Mignan et al. (2007), the contour limit of the stress shadow (i.e. the size of the stress shadow) is defined by the equation r⁎(t, σ = σ⁎) where σ⁎ approaches the failure stress σf and means that events are more likely to occur where r ≥ r⁎. The density per unit of space and time of events that compose the background seismicity is δ0b outside the stress shadow and δ−b inside the stress shadow, with δ0b N N δ−b.
44
A. Mignan / Tectonophysics 452 (2008) 42–50
Fig. 2. Schematic representation of the relation between the density of background seismicity δb and the stress σ. In an ideal case (left), events would only occur where the stress is positive (i.e. aftershocks) and regions of negative stress would remain inactive (i.e. quiescence). In a more realistic view (right), events are not influenced by very small stress changes (near-zero values) and occur randomly in space and time. These events compose the background seismicity (of density δ0b), which is generally observed in far field. Aftershocks (of density δ+b) occur where the stress is sufficient to trigger new events. In lobes of negative stress (i.e. stress shadow), there is a relative quiescence (of density δ−b) with only a few events being able to occur, such as δ+b N Nδ0b N N δ−b. The study of PAS only involves the spatiotemporal evolution of the stress limit σ⁎, which separates the background seismicity from erratic events in the stress shadow.
and time. These events compose the background seismicity of density δ0b, which is generally observed in far field. Aftershocks occur where the stress is sufficient to trigger new events. This triggering process increases the density of events to δ+b in lobes of positive stress. In lobes of negative stress (i.e. stress shadow), there is a relative quiescence with only a few events being able to occur for several reasons (e.g. fast stressing rate, transient effects, complex stress history). However, the density δ−b is supposed to be low with δ+b N N δ0b N N δ−b. Note that there is no reason for the density of events δ to evolve by step at given stress values. In particular for aftershocks, δ should increase continuously with increasing stress (I take the simple case where only static stress changes occur). However, this assumption has no consequences here since the study of PAS only involves the spatiotemporal evolution of the stress limit σ⁎, which separates the background seismicity (random in space) from erratic events in the stress shadow (seismicity noise that is also assumed random in space). I assume in this study that PAS (if it exists) is explained by the NC PAST. Other possible explanations (e.g. based on the theory of the Critical Point) are dismissed for the following reasons: (1) The NC PAST seems the most simple way to explain the powerlaw time-to-failure equation. It is based on simple and recognized concepts, such as elastic rebound and static stress transfer, and the origin of the power-law is purely geometrical. Therefore, I follow Ockham's Razor from which, of two equivalent theories or explanations, all other things being equal, the simpler one is to be preferred. (2) The NC PAST explains, in addition to PAS, other proposed seismicity precursors such as quiescence and the Mogi “doughnut” (e.g. Mogi, 1969). Note that this is already proposed in the Stress Accumulation model (King and Bowman, 2003). (3) Contrary to phenomena due to critical processes, which are difficult to catch in a complex dynamic system, PAS can be simply captured in space and time in the NC PAST view. This permits to determine how this precursor can be extracted, and thus to improve forecasting techniques. This aspect is used in the following simulations to test the power-law fit methodology (i.e. c-value). 3. Theoretical synthetic seismicity catalogues 3.1. Simulation methodology Theoretical synthetic catalogues are created based on the principles of the NC PAST. The background seismicity density δ0b and the noise ratio δ−b/δ0b are adjustable parameters that allow to test the consistency of PAS, depending respectively on the regional seismic
activity and on the effect of the stress shadow on seismicity. Coulomb stress calculations are realized using the source code AlmondX (courtesy of G. C. P. King), based on the dislocation theory in an elastic half-space (Okada, 1992). Events location x,y and time t are controlled by the stress contour limit σ⁎(x,y,t). The stress field σ(x,y,t) is calculated for a dip-slip fault (black segment) located in the center of a grid [0,1; 0,1]. The time interval is [t0 = 0; tf = 20]. For each time step ti (fixed to 1 year), events are random uniform in space and time (in this year interval) with a density δ0b outside the stress contour σ⁎(x,y,ti) and δ−b inside the stress contour σ⁎(x,y,ti) (Fig. 3, left panel). A more realistic representation of the behavior of background seismicity would take into account spatiotemporal clustering (i.e. ETAS model (e.g. Ogata and Zhuang, 2006) instead of a simple random uniform event distribution). However, since seismicity must be declustered in order to apply the power-law fit methodology, the latter approach seems more appropriate and straightforward. The stress field σ(x,y,ti) results from a displacement on the fault (imax − i)ΔD with ΔD being a constant displacement increment. At t0, an earthquake occurs on the fault with a displacement imaxΔD producing a stress shadow, which then decreases through time because of the component − iΔD (Fig. 3, right). This negative component corresponds to the linear superposition of a loading stress field by following the concept of back-slip model (Savage, 1983). The value of σ⁎ is fixed such that the last contour limit σ⁎(x,y,tf) is close to the fault segment (bL/10 with L the fault length), meaning that the fault is ready to fail (i.e. earthquake on the fault at t ~ tf, represented by a star on Fig. 3). Note that the value of σ⁎ depends on the fault mechanism, the displacement increment ΔD and the number of time steps ti. The relation r⁎(t, σ = σ⁎) that constrains the size of the stress shadow through time in one dimension (Mignan et al., 2007) can be applied to a two-dimensional case with r⁎(t, σ = σ⁎) valid along the gradient vectors of the stress field (two gradients represented on Fig. 3, right panel). However simulations are preferred since no simple analytical formulation can be given here. Fig. 4 shows an example of accelerating seismicity pattern (C = 0.29) observed in a theoretical synthetic catalogue. Events that compose the precursory signal are located in the region inside the stress contour σ⁎(t = t0), which corresponds to the maximum extent of the stress shadow, produced by the fault rupture at t = t0. The cumulative number of events λ increases as a power-law function through time, in agreement with Mignan et al. (2007). 3.2. Power-law fit methodology To detect accelerating seismicity patterns in a data set (here in theoretical synthetic catalogues), it is necessary to isolate the precursory events in a space/time window (also in a magnitude window in real cases, but not required here). Many studies use simple shapes such as circles (e.g. Brehm and Braile, 1998; Bowman et al., 1998; Robinson, 2000; Zoller et al., 2001) or ellipses (e.g. Papazachos et al., 2007) and a few other studies, based on the Stress Accumulation model (~NC PAST), use stress contours to define the region of interest (e.g. Mignan et al., 2006b; Bowman and King, 2001). In theoretical synthetic catalogues, PAS is located, by definition of the NC PAST, in the region delimited by a stress shadow. However in this study, since it is clearly impossible to determine the exact region of interest in real cases (i.e. stress contour depending on fault geometry and mechanism, on interactions with other active structures, etc…), PAS is identified in circular regions. This allows to obtain more realistic results from the optimization procedure. Fig. 5 (left panel) illustrates the optimization procedure employed for the systematic search of accelerating seismicity patterns in space and time. Patterns are searched for in concentric circles of different radii R centered on the fault and for different starting times tstart. The final time is fixed to tf. For each step (10 in space with 0.025 b R b 0.250
A. Mignan / Tectonophysics 452 (2008) 42–50
45
Fig. 3. Simulation of theoretical synthetic seismicity catalogues (NC PAST based). Events location x,y and time t are controlled by the stress contour limit σ⁎(x,y,t). The stress field σ(x,y,t) is calculated for a dip-slip fault (black segment) located in the center of a grid [0,1; 0,1]. The time interval is [t0 =0; tf = 20]. For each time step ti (fixed to 1 year), events are random uniform in space and time (in this year interval) with a density δ0b outside the stress contour σ⁎(x,y,ti) and δ−b inside the stress contour σ⁎(x,y,ti) (left panel). The value of σ⁎ is fixed such as the last contour limit σ⁎(x,y,tf) is close to the fault segment (bL/10 with L the fault length), meaning that the fault is ready to fail (i.e. earthquake on the fault at t ~tf, represented by a star). The relation r⁎(t, σ =σ⁎) that constrains the size of the stress shadow through time in one dimension (Mignan et al., 2007) can be applied to the two-dimensional case with r⁎(t, σ =σ⁎) being valid along the gradient vectors of the stress field (two gradients represented, right panel).
and 10 in time with 0 b tstart b 18), the c-value is calculated and the result plotted in the spatiotemporal search matrix or c-value matrix (right panel). In previous studies using the same optimization procedure (e.g. Sammis et al., 2004; Mignan et al., 2006a), the accelerating seismicity pattern is chosen as the signal with the lowest c-value (note that signals with a low c-value but composed of less than 10–15 events are generally dismissed). In theoretical synthetic catalogues, the optimal spatial extent and time duration of PAS are already known. The optimal starting time is t0 = 0 and in the example of Fig. 5, PAS is located in a circle of radius 0.10 b R b 0.15 (bold circles). R can only be approximated since no circle can fit the exact location of PAS (dark grey region). 4. Results
time in the interval 0 b t0 b10 (in theory, 0 b t0 b 3) and in a circle with 0.075 b R b 0.175 (in theory, 0.10 b R b 0.15). In addition, it must be noted that some instabilities can appear with a low c-value observed in more than one region (i.e. several blue patches). In such a case, it becomes difficult to determine the correct location of PAS. Note that Figs. A1–A5 clearly show the extreme variability of the c-value for different simulations with the same input parameters δ0b and δ−b/δ0b. The spatiotemporal variability/instability of the c-value can be explained by the fact that the power-law fit is very sensitive to the few events that occur in the first part of the acceleration (see Fig. 4). Because of the random behavior of events in space and time in the background, the starting point of the power-law curve is poorly controlled, leading to the variations observed among the c-value matrices. This tends to prove that the c-value does not permit to clearly determine the duration and spatial extent of PAS. This is due in part to the
4.1. Weakness of the power-law fit methodology Fig. 5 (right panel) shows an example of c-value matrix, which corresponds to a stack of 50 simulation results. In that particular case, the best acceleration (C b 0.3, blue region) is found for a starting time 0 b t0 b 3 and in a circle with a radius R comprised between 0.10 and 0.15 (bold circles on the left panel), which corresponds perfectly to the spatiotemporal window where the theoretical acceleration is located. This result is obtained thanks to the stacking method, which allows to increase the signal-to-noise ratio of the analysis, as proposed in some other forecasting techniques, such as Pattern Informatics (Holliday et al., 2005; Tiampo et al., 2002). Here the stacking method consists of creating a c-value matrix in which each value C(R, t0) corresponds to the mean of the values C(R, t0) computed from n simulations. The c-value matrices obtained for the 50 simulations (δ0b = 1000, δ−b/ δ0b = 0, see next section for different input parameters) are represented on Fig. A1 and four examples are represented on Fig. 6. Although PAS remains the same in theory (same input parameters δ0b and δ−b/δ0b), the random behavior of background seismicity that composes the precursory signal makes each pattern unique (example of Fig. 4). Differences do not change the shape of the acceleration at the first order but are sufficient to make significant changes on the power-fit result. On Fig. 6, the lowest c-value (C b 0.3) is obtained for a starting
Fig. 4. Example of accelerating seismicity pattern (C = 0.29) observed in a theoretical synthetic catalogue. Events that compose the precursory signal are located in the region inside the stress contour σ⁎(t = t0), which corresponds to the maximum extent of the stress shadow, produced by the fault rupture at t = t0. The cumulative number of events λ increases as a power-law function through time, in agreement with Mignan et al. (2007). The power-law fit and linear fit used to calculate the c-value are represented.
46
A. Mignan / Tectonophysics 452 (2008) 42–50
Fig. 5. Method employed for a systematic search of accelerating seismicity patterns in space and time. Left panel: Patterns are searched in circles of different radii centered on the fault and for different starting times. The final time is fixed to tf. Right panel: Spatiotemporal search matrix or c-value matrix. C-values are computed for 10 starting times in the interval [0; 18] and for 10 radii in the interval [0.025; 0.250]. In this particular example, the result corresponds to a stack of c-value matrices obtained from 50 simulations (δ0b = 1000, δ−b/δ0b =0.).
inherent weakness of the power-law fit but also to the random behavior of events that compose PAS. If the first issue (linked to the fit methodology) could only be resolved by adopting a new approach to extract PAS (not the purpose of this article), the second issue (linked to the random behavior of events) can be resolved by stacking many data sets. 4.2. Consistency of PAS As already discussed before, the existence of AMR (~PAS) is usually questioned on the basis that the duration and spatial extent of AMR (~PAS)
are highly variable and do not seem to be linked to any characteristic of the subsequent mainshock. This argument seems at present irrelevant since this variability is also observed in theoretical synthetic catalogues where PAS is known to exist. The other argument against the existence of AMR (~PAS) is that it is not observed systematically before large earthquakes. This argument cannot be refuted by the results of the previous simulations (δ0b =1000, δ−b/δ0b =0.), where a clear acceleration (i.e. low c-value) is found in the majority of cases (see Figs. 6 and A1). Therefore some new simulations are realized with different input parameters (δ0b and δ−b/δ0b) in order to test the consistency of PAS. First,
Fig. 6. Four examples of c-value matrix. Although the input parameters of the different simulations are the same (i.e. δ0b = 1000, δ−b/δ0b = 0), the power-law fit result is highly variable. This is due in part to the inherent weakness of the power-law fit methodology but also to the random behavior of background seismicity that composes PAS. See Fig. A1 for the c-value matrices of 50 simulations.
A. Mignan / Tectonophysics 452 (2008) 42–50 Fig. 7. Role of the background seismicity density δ0b in the formation of PAS (noise ratio δ−b/δ0b fixed to zero). From left to right, δ0b = 1000, δ0b = 500 and δ0b = 100 (i.e. number of events in the unit square [0,1; 0,1], per year). For each value of δ0b, the upper panel shows an example of theoretical PAS. The central panel shows the probability density function (PDF) of the corresponding c-value, for 100 simulations. The lower panel shows the c-value matrix obtained after stacking the results of 50 simulations (All c-value matrices are represented in Figs. A1–A3 respectively). 47
48 A. Mignan / Tectonophysics 452 (2008) 42–50 Fig. 8. Role of the noise ratio δ−b/δ0b in the formation of PAS (background density δ0b fixed to 1000). From left to right, δ−b/δ0b = 0.0, δ−b/δ0b = 0.2 and δ−b/δ0b = 0.5. For each value of δ−b/δ0b, the upper panel shows an example of theoretical PAS. The central panel shows the probability density function (PDF) of the corresponding c-value, for 100 simulations. The lower panel shows the c-value matrix obtained after stacking the results of 50 simulations (All c-value matrices are represented in Figs. A1, A4 and A5 respectively).
A. Mignan / Tectonophysics 452 (2008) 42–50
the possible role of a decrease of the regional seismic activity on PAS is examined by changing the value of δ0b (noise ratio δ−b/δ0b fixed to zero). Fig. 7 illustrates the evolution of PAS for δ0b = 1000, 500 and 100 (i.e. number of events in the unit square [0,1; 0,1], per year). For each value of δ0b, the upper panel shows an example of theoretical PAS, which is defined from its exact spatiotemporal extent. The number of events in the signal is λf. The central panel shows the probability density function (PDF) of the corresponding c-value, for 100 simulations. The lower panel shows the c-value matrix obtained after stacking the results of 50 simulations (all c-value matrices are represented in Figs. A1–A3 respectively). It can be observed that a change of the background seismicity density δ0b has a clear impact on the consistency of PAS. The lower the density δ0b, the more the precursory signal tends to disappear. On the c-value PDF (central panel), the coefficient of variation increases and the mean value shifts to higher c-values. This phenomenon is not a consequence of the weakness of the power-law fit since it can still be observed after stacking many c-value matrices (lower panel) but is due to the way PAS forms in the NC PAST. In that view, an ideal acceleration would form if δ0b tends to infinity. Mignan et al. (2007) defined the power-law time-to-failure equation for a continuous quantity, linked to the continuous decrease of the size of a stress shadow. Since the precursory signal is composed of punctual events, some information is irremediably lost. With δ0b decreasing, the signal fluctuates more and more around the ideal power-law curve, as observed on the upper panel of Fig. 7. Fig. 8 shows the impact of an increase of the noise ratio δ−b/δ0b on PAS for δ−b/δ0b = 0.0, 0.2 and 0.5 (with δ0b fixed to 1000). In the same way, the higher the noise ratio δ−b/δ0b, the more the precursory signal tends to disappear. In that case, the theoretical PAS is still present (clear for δ0b = 1000) but a linear trend is superimposed to the power-law. Therefore the power-law fit methodology is less efficient. With δ−b/ δ0b = 0.2, the best c-value obtained in the spatiotemporal search matrix (lower panel) increases from 0.3 to 0.5 and with δ−b/δ0b = 0.5, the signal becomes very weak (C ~ 0.7), which would not be considered as a reliable signal in a real case. These two tests permit to explain why PAS (or AMR) is not systematically observed before large earthquakes. If the NC PAST is correct, PAS would only be clearly observed in regions of high seismic activity. Moreover, quiescence phenomena are not systematically observed in stress shadows. Hardebeck et al. (1998) showed that only 65% of aftershocks occurred in regions of stress increase related to the 1992 Landers earthquake, 60% to the 1994 Northridge earthquake. Therefore between 35 and 40% of all aftershocks occurred in regions of decreased stress, part of these aftershocks certainly occurred in the stress shadow delimited by the negative stress value σ⁎. This tends to prove that the noise ratio δ−b/δ0b might be high in real seismicity catalogues, leading to precursory patterns difficult to catch. 5. Discussion The purpose of this article has not been to prove or to verify that PAS exists but to show, assuming that the NC PAST is correct, that one can explain why PAS is not observed systematically before large earthquakes and why the spatial extent and duration of the signal seem unstable using the c-value. It is important to note that Bowman et al. (1998) clearly formalized the AMR concept into a testable hypothesis thanks to the c-value (noted by Hardebeck et al., submitted for publication). The c-value provides a simple, computationally inexpensive parameter to search for AMR/PAS patterns and thus has been used in many studies since Bowman et al. (1998) (see review in Mignan, 2007). Nowadays, authors rarely make the distinction between c-value and AMR, leading irremediably to some misconceptions. I have shown in this study that an unstable c-value does not mean that AMR/PAS does not exist. I have demonstrated why the signal can
49
be weak and what parameters control its appearance (i.e. the background seismic activity δ0b and the background seismicity noise ratio δ−b/δ0b). These results lead to three important conclusions: (1) The c-value is inadequate for robust systematic prospective forecasts. The c-value only gives consistent results for regions of high background seismic activity (δ0b high) where stress shadows are clearly defined (δ−b/δ0b low). This case is rarely encountered in real seismicity catalogues. It is interesting to note that stacking c-value matrices permits to improve the power-law fit results. However this can only be done in retrospect and requires that PAS patterns due to different large earthquakes have the same spatiotemporal extent, which is also unlikely. (2) PAS/AMR possibly exists if the NC PAST is correct. This new theory makes specific predictions for the spatiotemporal extent of PAS and explains why the signal can appear or not. Future studies based on the NC PAST should focus on some new methodologies that would permit to extract PAS in a better way. The c-value has been proved to be simple and useful in the past but new results show that more complicated algorithms are necessary to possibly forecast large earthquakes. (3) Another alternative would be to improve the seismic network detection threshold. The NC PAST shows indeed that it is the number of events which is important and not the magnitude for the appearance of PAS. This could be done by increasing the density of simple one component narrow band seismometers in shallow boreholes. At present the trend in networks is to have fewer 3 component instruments with the result that detection threshold has declined in many places such as South California and Greece (G. C. P. King, pers. comm.). Acknowledgements The author acknowledges helpful reviews by Kristy Tiampo and Geoffrey King as well as Editor Tom Parsons. The author also thanks Geoffrey King for the source code AlmondX used for the Coulomb stress calculations. Appendix A. Supplementary data Supplementary data associated with this article can be found, in the online version, at doi:10.1016/j.tecto.2008.02.010. References Bowman, D.D., King, G.C.P., 2001. Accelerating seismicity and stress accumulation before large earthquakes. Geophys. Res. Lett. 28, 4039–4042. doi:10.1029/2001GL013022.2001. Bowman, D.D., Ouillon, G., Sammis, C.G., Sornette, A., Sornette, D., 1998. An observational test of the critical earthquake concept. J. Geophys. Res. 103, 24,359–24,372. Brehm, D.J., Braile, L.W., 1998. Intermediate-term earthquake prediction using precursory events in the New Madrid Seismic Zone. Bull. Seismol. Soc. Am. 88, 564–580. Bufe, C.G., Varnes, D.J., 1993. Predictive modeling of the seismic cycle of the greater San Francisco Bat region. J. Geophys. Res. 98, 9871–9883. Hardebeck, J.L., Nazareth, J.J., Hauksson, E., 1998. The static stress change triggering model: constraints from two southern California aftershock sequences. J. Geophys. Res. 103, 24,427–24,437. Hardebeck, J.L., Felzer, K.R., Michael, A.J. (submitted for publication), Rigorous observational tests contradict the accelerating moment release hypothesis. Holliday, J.R., Nanjo, K.Z., Tiampo, K.F., Rundle, J.B., Turcotte, D.L., 2005. Earthquake forecasting and its verification. Nonlinear Process Geophys. 12, 965–977. Jaume, S.C., Sykes, L.R.,1999. Evolving towards a critical point: a review of accelerating seismic moment/energy release prior to large and great earthquakes. Pure Appl. Geophys. 155, 279–306. King, G.C.P., Bowman, D.D., 2003. The evolution of regional seismicity between large earthquakes. J. Geophys. Res. 108. doi:10.1029/2001JB000783. Mignan, A., 2007. The Stress Accumulation model: accelerating moment release & seismic hazard. Adv. Geophys. 49. doi:10.1016/S0065-2687(07)49002-1. Mignan, A., King, G.C.P., Bowman, D.D., Lacassin, R., Dmowska, R., 2006a. Seismic activity in the Sumatra-Java region prior to the December 26, 2004 (Mw = 9.0–9.3) and March 28, 2005 (Mw = 8.7) earthquakes. Earth Planet. Sci. Lett. 244, 639–654.
50
A. Mignan / Tectonophysics 452 (2008) 42–50
Mignan, A., Bowman, D.D., King, G.C.P., 2006b. An observational test of the origin of accelerating moment release before large earthquakes. J. Geophys. Res. 111. doi:10.1029/2006JB004374. Mignan, A., King, G.C.P., Bowman, D.D., 2007. A mathematical formulation of accelerating moment release based on the stress accumulation model. J. Geophys. Res. 112. doi:10.1029/2006JB004671. Mogi, K., 1969. Some features of recent seismic activity in and near Japan (2), Activity before and after large earthquakes. Bull. Earthq. Res. Inst. Univ. Tokyo 47, 395–417. Ogata, Y., Zhuang, J., 2006. Space–time ETAS models and an improved extension. Tectonophysics 413, 13–23. Okada, Y., 1992. Internal deformation due to shear and tensile faults in a half-space. Bull. Seismol. Soc. Am. 82, 1018–1040. Papazachos, B.C., Karakaisis, G.F., Papazachos, C.B., Scordilis, E.M., 2007. Evaluation of the results for an intermediate-term prediction of the 8 January 2006 Mw = 6.9 Cythera earthquake in the Southwestern Aegean. Bull. Seismol. Soc. Am. 97, 347–352. doi:10.1785/0120060075. Reid, H.F., 1910. The mechanics of the earthquake, the California earthquake of April 18, 1906. Report State Investig. Comm., vol. 2. Carnegie Inst., Washington.
Robinson, R., 2000. A test of the precursory Accelerating Moment Release model on some recent New Zealand earthquakes. Geophys. J. Int. 140, 568–576. Sammis, C.G., Sornette, D., 2002. Positive feedback, memory, and the predictability of earthquakes. Colloq. Proc. - Natl. Acad. Sci. 99, 2501–2508. Sammis, C.G., Bowman, D.D., King, G.C.P., 2004. Anomalous seismicity and accelerating moment release preceding the 2001 and 2002 earthquakes in northern Baja California, Mexico. Pure Appl. Geophys. 161, 2369–2378. Savage, J.C., 1983. A dislocation model of strain accumulation and release at a subduction zone. J. Geophys. Res. 88, 4984–4996. Tiampo, K.F., Anghel, M., 2006. Introduction to special issue: critical point theory and space–time pattern formation in precursory seismicity. Tectonophysics 413, 1–3. Tiampo, K.F., Rundle, J.B., McGinnis, S., Klein, W., 2002. Pattern dynamics and forecast methods in seismically active regions. Pure Appl. Geophys. 159, 2429–2467. Zoller, G., Hainzl, S., Kurths, J., 2001. Observation of growing correlation length as an indicator for critical point behavior prior to large earthquakes. J. Geophys. Res. 106, 2167–2176.