CC BY
17
RUSSIAN JOURNAL OF EARTH SCIENCES, VOL. 19, ES6003, doi:10.2205/2019ES000695, 2019
Epidemic type aftershock sequence exponential productivity
S. V. Baranov1, A. D. Gvishiani2'3, C. Narteau4, and P. N. Shebalin5
Received 20 September 2019; accepted 22 November 2019; published 30 November 2019.
Studying the hierarchical structure of the aftershock sequence of the three largest earthquakes of the last decade, we show that the number of offspring events counted in a fixed magnitude band with respect to the magnitude of the parent events follows an exponential distribution. Such an exponential productivity law is coherent with the exponential decays inferred from largest earthquakes worldwide. Epidemic Type Aftershock Sequences (ETAS) are the most popular stochastic models of seismicity and they are all based on a Poisson distribution of the earthquake productivity with a pronounced non-zero mode. We construct here an alternative model incorporating an exponential productivity law. For the three aftershock sequences. we estimate parameters of both models using aftershocks occuring during the first 2 days. We simulate a set of synthetic earthquake catalogues for both models and compare the average cumulated number of events with respect to time. In all cases, the ETAS model overestimates the number of events in the interval from 2 to 365 days. For the same time period, the exponential ETAS model gives a satisfactory cumulative number of events. We conclude that exponential distribution of the earthquake productivity seems to be an important property of the seismic relaxation process. KEYWORDS: Strong earthquakes; aftershocks; productivity; ETAS model; EP model; forecast.
Citation: Baranov, S. V., A. D. Gvishiani, C. Narteau, and P. N. Shebalin (2019), Epidemic type aftershock sequence exponential productivity, Russ. J. Earth. Sci, 19, ES6003, doi:10.2205/2019ES000695.
Introduction the way in which a sequence of earthquakes develops over space and time, we may consider evAn important feature of seismicity is the oc- ery event as the trigger for the perturbation of the currence of space-time clusters demonstrating that state of stress in some area around the event. The earthquakes interact with each other. Focusing on zones of such a perturbation intersect in time and
__space, [Gvishiani et al., 2013a; 2013b; 2016], and
lKola Branch Federal Research Center Ge°physical each new event may be considered as an "offspring"
Survey, Russian Academy of Sciences, Apatity, Russia of an preceding earthquakes. For a Poisson pro-
2 Geophysical Center of Russian Academy of Sciences, , , , ,
Moscow Russia cess such epidemic behavior may be described as
3Institute of the Physics of the Earth, Russian a superposition of the processes initiated by every
Academy of Sciences, Moscow, Russia event. Since the advent of epidemic models of seis-
4Universite de Paris, Institut de physique du globe micity, the productivity has become a major issue
de Paris, CNRS, F-75005 Paris, France. because it determines the increase in seismicity rate
Insti^te of Earthquake ^dictton Theory and after each earthquake [Helmstetter, 2002; Kagan, Mathematical Geophysics, Moscow, Russia
1981; Ogata,, 1989]. In these models the number
„ . ,,„„-,„, „ n i • i ^ ^ n ao of events triggered by a magnitude m earthquake
Copyright 2019 by the Geophysical Center RAS. . , -rf . r
is considered to vary as a Poisson process of rate http://rjes.wdcb.ru/doi/2019ES000695-res.html J
ES6003
1 of 8
(N(m)) = K 10am. The value a was estimated in a range from 0.5 to 2 [Console et al., 2003; Hainzl and Marsan, 2008; Hainzl et al., 2013; Wang et al., 2010; Werner and Sornette, 2008; Zhuang et al., 2004; Zhuang et al., 2005]. This value is usually close to the observed 6-value, the slope of the earthquake-size distribution [Helmstetter, 2003].
However, these estimates remain uncertain due to the difficulty of isolating the relative contributions of successive events in a sequence. This task remains in early stage despite the diversity of declustering methods implemented in the past.
Two stochastic approaches have been suggested to find causal links within cascades of triggered seismicity. A first approach is to separate the branching structure of earthquake sequences from the background rate using an iterative algorithm based on maximum likelihood estimation of the parameters of an epidemic model of seismicity ETAS [Zhuang et al. 2002]. A second approach is model independent.
A linear contributions of each earthquake to the overall seismicity rate is assumed [Marsan and Helmstetter. 2008]. Those two approaches suppose that every new earthquakes is an "aftershock" of all preceding events, and the goal is to estimate the impact of each preceding event on each subsequent event in terms of probability. An alternative approach [Zaliapin and Ben-Zion, 2013; Zaliapin et al., 2008] goes directly to consider a tree of events in which each event may be a "parent" of several later events, but it can be an "offspring" of only one earlier event. Technically the "parent" is found as a "nearest neighbor" using proximity functions in time-space-magnitude domains [Baiesi and Paczuski, 2004; Zaliapin et al., 2008]. Here we shall call "parent" events triggering events, and their "offsprings" triggered events.
All these methods confirm the dependency of the productivity on the magnitude of the triggering event. However, it was found that within this dependency there is a huge variability in the number of triggering events in the seismic catalogs [Marsan and Helmstetter, 2017]. Recently it was found that unexpectedly this variability may be described by exponential distribution [Shebalin et al. 2018] with maximum at 0. This result was obtained for the global statistics of the number of aftershocks from earthquakes of magnitude 6.5 and above. Aftershocks of magnitude above a threshold, defined as
the magnitude of the main shock minus 2, were counted.
Here we study a distribution of the productivity in a tree of aftershocks for three largest earthquakes of the last decade: 11 March 2011, Mw=9.1, Tohoku earthquake; 27 February 2010, Mw=8.8, Chile earthquake; 11 April 2012, Mw=8.6, Sumatra earthquake. We find that for each sequence the distribution of the productivity also tends to an exponential form. Thus, exponential distribution of earthquake productivity seems to be a general property of seismicity. Exponential distribution of the real number of triggered events contradicts to usually expected form of Poisson distribution. To better understand this issue we compare two models for the three aftershock sequences.
Earthquake Productivity
We define earthquake productivity using "delta-analysis" [Zaliapin and Ben-Zion, 2013; Zaliapin et al., 2008], in a magnitude band of a given width AM relative to the magnitude of each triggering event. Thus, the productivity is a property of each earthquake.
Tree of Clustered Earthquakes
To count the productivity values we decompose the earthquake catalogue into a hierarchical tree of pairwise links triggering-triggered events. For each pair of earthquakes {i, j}, we compute the proximity function[Baiesi and Paczuski, 2004],
=
{
tij (rZJ )df 10
—bm,i
for tij > 0, for tij < 0.
(1)
where tij = tj — t
is the inter-event time, rij the spatial distance between the epicenters, mi the magnitude of event i, df the fractal dimension of the epicenter distribution and b the slope of the earthquake-size distribution.
Using the proximity function (1) for each event we find the preceding nearest-neighbor. In case the ^ value exceeds a threshold the event is considered as a background event because it has no triggering event. For each sequence we computed rj0 using the original technique of [Zaliapin
and Ben-Zion, 2013; Zaliapin et al., 2008], approximating the distribution of the nearest-neighbor values log(rj) by a mixture model of two Gaussian distributions, one modeling independent events, the other causally-related events. We select the value for which the two types of errors compensate each other: same probability of having causally-related events with ^ > and independent events with ^ <rq0.
Accordingly, all clusters are built from a primary triggering event. There is a single path from any earthquake in a cluster to the corresponding primary event. Primary events by definition are "background events". Background events with the largest magnitude in the cluster are mainshocks. But a triggered event may have a larger magnitude than its triggering events. In this case the main shock is not a primary event and nor background event. Those definitions describe an aftershock sequence as a hierarchical tree of events. This is slightly different from a standard definition of the foreshocks - main shock - aftershock series in which the majority of the aftershocks (and also foreshocks) are linked directly to the main shock. In order to avoid the excessive influence of the definition of the proximity function on our analysis, for each of the three series we consider not a single hierarchical tree of events, including the mentioned earthquakes, but all earthquakes in a simple spatial region with a form of stadium [Baranov and Shebalin, 2017] in a time interval of 365 days after the earthquake.
For each triggering event, we count the number N of triggered events at the lower hierarchical level using a relative magnitude threshold AM (i.e. ^triggering — ^triggered < AM). This number N is defined as the productivity.
Distribution of the Productivity
ing earthquakes nor on the width of the considered magnitude band. For the exponential distribution the clustering factor is an important parameter, as it is a single parameter of this distribution.
Here we concentrate on the productivity distribution within aftershock sequences of large earthquakes. It is important to know whether the exponential form is also characteristic in much more homogeneous conditions in comparison to the global variability of the aftershock sequences considered earlier [Shebalin et al. 2018].
The estimated completeness magnitude for the Sumatra and Chile sequences is 4.5, for Tohoku 5.0. For the analysis we use the value AM = 1.5. This choice allowed to choose the minimum magnitude of triggering earthquakes Mtr 6.0 for the Sumatra and Chile sequences and 6.5 for the To-hoku sequence. Figure 1 shows the histograms of the productivity in comparison to exponential and Poisson distribution. We note that exponential distribution is the distribution of a real value, and the actual productivity is an integer. We may interpret this challenge by supposing that the productivity is an internal property of each earthquake, similar to its magnitude. Specific realizations of the productivity is an integer. It is natural to suppose that the specific values have Poisson distribution with a rate equal to the "internal" productivity. The existing epidemic models of seismicity imply that the internal productivity is constant. The difference between the "internal" and the real productivity may explain some distortion of the exponential distribution. However it is clear from the figure that for all the three sequences exponential distribution is preferable in comparison to Poisson distribution.
Two Alternative Epidemic Models of Seimicity
The distribution of the number of triggered events We test here whether the found property is im-for an earthquake population is defined as the pro- portant for modeling the seismicity. We compare ductivity distribution with a mean denoted Ло(AM), two epidemic models, the classic ETAS model and we call the clustering factor. The productivity may a similar model in which the constant internal pro-vary from place to place, and also from sequence ductivity is replaced by a random internal producto sequence. Recently it was found [Shebalin et tivity with exponential distribution. Using an in-al. 2018] that in a global scale the productivity terval of 2 days right after the large earthquakes. has a distribution of the exponential shape with Using those estimates we construct two versions maximum at 0. It was also shown that this shape of a synthetic catalog. Repeating simulations many does not depend on the magnitude of the trigger- times, we calculate an average cumulative number
Figure 1. Earthquake productivity for three aftershock sequences. (a) - 27 February 2010, Mw=8.8, Chile earthquake, (b) - 11 March 2011, Mw=9.1, Tohoku earthquake, (c) - 11 April 2012, Mw=8.6, Sumatra earthquake. Dots show the number of triggered events for M > 6.5 (Tohoku) and M > 6.0 (Chile and Sumatra) triggering events using a relative magnitude threshold AM = 1.5. The solid line is the exponential law with parameter Ao, the clustering factor. The histogram shows the Poisson distribution with parameter A0.
of events as a function of time in the interval 0 to 365 days.
Finally we compare the results for the two models with the real data.
ETAS Model
The classic ETAS model [Kagan and Knopoff, 1981; Ogata, 1989] considers seismicity as a non-stationary Poisson process with a rate described by the equation (2).
10«(MS-Mc) A ETAS (t)= » + K^(t _ ti + C)P,
(2)
where ti and Mi are the time and the magnitude of the f-th earthquake, K0, a, and p are parameters, M0 magnitude threshold for counting events. Parameters and p describe temporal decay of the triggered events according to the empirical Omori-Utsu law [Utsu, 1965; Utsu, 1970]. Parameter a, as discussed above, is usually close to the 6-value of the Gutenberg-Richter relation. Parameter ^ is the background seismicity rate, which is assumed constant.
To estimate parameters of the model we use a standard maximum likelihood procedure [Ogata, 1989]. Earthquake catalogs are not complete right after large earthquakes. To minimize the impact of this effect we omit in the analysis the interval upto 0.05 days after the earthquake.
For the synthetic catalog simulations there are two equivalent ways [Zhuang et al., 2004]. The first is to treat the ETAS model as a single point process, with the probability of an event at each point in time reflected in the conditional intensity which contains a component of background and a component of triggered seismicity contributed by all past events in the history of the process. The second way is to simulate the background events as a Poisson process, and then recursively simulate the aftershocks resulting from each of these background events in turn. Finally, all events are combined and put into the correct temporal order of occurrence. The first way is not appropriate for our modification of the ETAS model, because the internal productivity is random. For this reason we apply here the second way. In both methods the magnitude of events is simulated independently as a random number above M0 with distribution defined by Gutenberg-Richter relation.
Table 1. ETAS and EP models parameters estimated using data for 2 days after mainshock
Mainshock 6-value Model parameters
Chili, 2010, M8.8 b = 1.36 ETAS: ß = 0.65, c = 0.08, p = 1.26, KQ = 0.028, a = 1.82
EP: AM = 1.5, ßEp = 0.1, A = 0.2, c = 0.034, p = 1.01 Tohoku, 2011, M9.1 b = 1.06 ETAS: ß = 0.81, c = 0.17, p = 1.56, K0 = 0.04, a = 1.74
EP: AM = 1.5, ßEP = 0.2, A = 1, c = 0.04, p = 1.02 Sumatra, 2012, M8.6 b = 1.26 ETAS: ß = 0.43, c = 0.1, p = 1.52, K0 = 0.01, a = 1.846
EP: AM = 1.5, ßEP = 0.1, A = 0.7, c = 0.05, p = 1.36
ETAS Model with Random Productivity with Exponential Distribution
We suggest here a simple modification of the ETAS model we shall call EP model. This model modifies the recursive definition of the ETAS model described above. The model represents seismicity as a sum of Poisson processes with the decay according to the Omori law, initiated by background events. Each new event also initiates a similar decaying process. For the synthetic catalogue simulations for each new event we randomly generate its internal productivity Aj according to the exponential distribution
1
p(Ai) = — eA*/A° A0
Then we simulate each branch of the non-stationary process with the rate defined by the equation
(3):
AEPi (t) =
A
(3)
(1 + (i — U)/C)P'
In simulations we use Bayesian estimates of the parameters c and p with Gaussian priors [Shebalin and Baranov, 2019] in the interval (0.05,2) days after the large earthquakes. The clustering factor A0 is also estimated in this interval simply as the average productivity.
This number is corrected to the interval of 365 days using a multiplier
U =
U(0.05, 365) U(0.05, 2) :
where U(a,b) = £ (1 + x/c)-pdx. Like in ETAS model, magnitudes of events are independently assigned according to the Gutenberg-Richter relation. We use a Bayesian estimate of the 6-value with Gaussian priors [Shebalin and Baranov, 2019].
Data
We used data from ANSS ComCat earthquake catalog provided by USGS. Aftershocks of M8.8 Chile earthquake of 2010 were taken for a year after the mainshock from the circle of radius 900 km surrounding its epicenter. Aftershocks of M9.1 Tohoku earthquake of 2011 were taken for a year after the mainshock from the circle of radius 1000 km surrounding its epicenter. Aftershocks of M8.6 Sumatra were taken for a year after the mainshock from the circle of radius 700 km surrounding its epicenter. We did not apply any other special aftershock selection procedure.
Results and Discussion
For each of the three sequences we have performed 2500 simulations of the synthetic catalogs using the two models. We estimated ETAS and EP models parameters using data for 2 days after the mainshocks (Table 1). The average cumulative number of event has been calculated as a function of time. Results are shown in Figure Figure 2.
We see that the ETAS model gives a drastic over-estimation of the earthquake rates at later times, while the EP model is quite acceptable in the whole forecasting period from 2 to 365 days. Two major reasons explain this as we can suppose. First, values following the exponential distribution are statistically smaller than values of the Poisson distribution with the same mean value. Thus, the EP model should predict smaller rates in comparison to the ETAS model. Second, in all three cases the background rate of the ETAS model estimated at the beginning of the sequence is much higher than estimated in the interval of 365 days. In the EP model the calculations of the background rate give
ES6003 BARANQY ET AL.: EPIDEMIC TYPE AFTERSHOCK SEQUENCE ES6003
Figure 2. Cumulative number of seismic events for three aftershock sequence, real and forecasted by models ETAS and EP. (a) - 27 February 2010, Mw=8.8, Chile earthquake, (b) - 11 March 2011, Mw=9.1, Tohoku earthquake, (c) - 11 April 2012, Mw=8.6, Sumatra earthquake. Circles show the real cumulative number of event of magnitude M > 6.5 (Tohoku) and M > 6.0 (Chile and Sumatra), dashed line the average for 2500 catalogue simulation of the ETAS model, solid line average for 2500 catalogue simulation of the EP model. Parameters for both models have been estimated in the interval (0.05,2) days afters the corresponding large eartquakes, their value are given in Table 1.
ES6003
BARANOY ET AL.: EPIDEMIC TYPE AFTERSHOCK SEQUENCE
ES6003
similar values for the beginning and for the whole sequence.
One of interesting results of this paper is exponential form of the productivity distribution within single aftershock sequences from large earthquakes. Earlier this property was established in global and regional scales [Shebalin et al. 2018; Shebalin and Baranov, 2019].
The obtained results demonstrate advantages of the suggested modification of the ETAS model. It is too early to come to final conclusions about those two model, however we obviously may conclude that exponential shape of the productivity distribution is an important property of the seismic process.
Acknowledgments. AG acknowledges financial support from the Presidium of the Russian Academy of Sciences Basic Research Program No. 2, project "Complex physical and mathematical modelling of circumpolar ionosphere and the discrete mathematics techniques in the tasks of digital analysis of geophysical data using supercomputers". SB developed software for the data analysis with support of Russian Foundation of Basic Researches, project 19-05-00812.
References
Baiesi, M., M. Paczuski (2004), Scale-free networks of earthquakes and aftershocks, Phys. Rev. E, 69, 066106, Crossref Baranov, S. V., V. A. Pavlenko, P. N. Shebalin (2019), Forecasting Aftershock Activity: 4.
Estimating the Maximum Magnitude of Future Aftershocks, Izv., Phys. Solid Earth, 55, 548-562, Crossref
Baranov, S. V., P. N. Shebalin (2017), Forecasting aftershock activity: 2. Estimating the area prone to strong aftershocks, Izv., Phys. Solid Earth, 53, 366384, Crossref Console, R., M. Murru, A. M. Lombardi (2003), Refining earthquake clustering models, J. Geophys. Res.: Solid Earth, 108, Crossref Gvishiani, A., M. Dobrovolsky, S. Agayan, et al. (2013a), Fuzzy-based clustering of epicenters and strong earthquake-prone areas, Environmental Engineering and Management Journal, 12, No. 1, 1-10, Crossref
Gvishiani, A., B. Dzeboev, S. Agayan (2013b), A new approach to recognition of the earthquake-prone areas in the Caucasus, Izvestiya, Physics of the Solid Earth, 49, No. 6, 747-766, Crossref Gvishiani, A., B. Dzeboev, S. Agayan (2016),
FCAZm intelligent recognition system for locating areas prone to strong earthquakes in the Andean and Caucasian mountain belts, Izvestiya.. Physics of the Solid Earth, 52, No. 4, 461-491, Crossref Hainzl, S., D. Marsan (2008), Dependence of the omori-utsu law parameters on main shock magnitude: Observations and modeling, J. Geophys. Res.: Solid Earth, 113, Crossref Hainzl, S., O. Zakharova, D. Marsan (2013), Impact of aseismic transients on the estimation of aftershock productivity parameters, Bull. Seismol. Soc. Am, 103, 1723-1732, Crossref Helmstetter, A. (2003), Is earthquake triggering driven by small earthquakes?, Phys. Rev. Lett., 91, 058501, Crossref Helmstetter, A., D. Sornette (2002), Subcritical and supercritical regimes in epidemic models of earthquake aftershocks, J. Geophys. Res.: Solid Earth, 107, 2237, Crossref Holschneider, M., G. Zoller, S. Hainzl (2011), Estimation of the maximum possible magnitude in the framework of the doubly-truncated gutenberg-richter model, Bull. Seimol. Soc. Am, 112, 1649-1659, Crossref
Kagan, Y. Y., L. Knopoff (1981), Stochastic
synthesis of earthquake catalogs, J. Geophys. Res.: Solid Earth, 86, 2853-2862, Crossref Marsan, D., A. Helmstetter (2017),. Single-link cluster analysis as a method to evaluate spatial and temporal properties of earthquake catalogues, J. Geophys. Res.: Solid Earth, 122, 5544-5560, Crossref Marsan, D., O. Lengline (2008), Extending
earthquakes' reach through cascading, Science, 319, 1076-1079, Crossref Ogata, Y. (1989), Statistical model for standard seis-micity and detection of anomalies by residual analysis, Tectonophysics, 169, 159-174, Crossref Shebalin, P. N., S. V. Baranov (2019), Forecasting Aftershock Activity: 5. Estimating the Duration of a Hazardous Period, Izv., Phys. Solid Earth, 55, 719-732, Crossref Shebalin, P. N., S. V. Baranov, B. A. Dzeboev (2018), The Law of the Repeatability of the Number of Aftershocks, Dokl. Earth Sc., 481, 963-966, Crossref Utsu, T. (1965), A method for determining the value of b in a formula logn = a _ bM showing the magnitude-frequency relation for earthquakes (with English summary), Geophys Bull. Hokkaido Univ., 13, 99-103.
Utsu, T. (1970), Aftershocks and earthquake statistics (ii): Further investigation of aftershocks and other earthquake sequences based on a new classification of earthquake sequences, J. Faculty Sci., Hokkaido University, Ser. VII (Geophys.), 3, 197-266. Wang, Q., F. P. Schoenberg, D. D. Jackson (2010), Standard errors of parameter estimates in the etas model, Bull. Seismol. Soc. Am., 100, 1989-2001, Crossref
Werner, M. J., D. Sornette (2008), Magnitude uncertainties impact seismic rate estimates, forecasts,
and predictability experiments, J. Geophys. Res.: Solid Earth, 113, Crossref Zaliapin, I., Y. Ben-Zion (2013), Earthquake clusters in southern california i: Identification and stability, J. Geophys. Res.: Solid Earth, 118, 28472864, Crossref Zaliapin, I., A. Gabrielov, V. Keilis-Borok, H. Wong (2008), Clustering analysis of seismicity and aftershock identification, Phys. Rev. Lett., 101, 018501, Crossref
Zhuang, J., C.-P. Chang, Y. Ogata, et al. (2005), A study on the background and clustering seismicity in the taiwan region by using point process models, J. Geophys. Res.: Solid Earth, 110, Crossref
Zhuang, J., Y. Ogata, D. Vere-Jones (2002), Stochastic declustering of space-time earthquake occurrences, J. Am. Stat. Ass., 97, 369-380, Crossref Zhuang, J., Y. Ogata, D. Vere-Jones (2004), Analyzing earthquake clustering features by using stochastic reconstruction, J. Geophys. Res.: Solid Earth, 109, Crossref
Corresponding author:
P. N. Shebalin, Institute of Earthquake Prediction Theory and Mathematical Geophysics, Moscow, Russia. (p.n.shebalin@gmail.com)
