Application of Time-to-Failure Model to Precursory Seismic
Activation on Kamchatka and in Italy
IUGG 99 Report (IUGG 99 Abstracts, p. A.141)
R. Di Giovambattista
Istituto Nazionale di Geofisica, Via di Vigna Murata, 605, I-00143 Roma, Italy. E-mail: digiovam@marte.ingrm.it
G. A. Sobolev
Institute of the Physics of the Earth, B.Gruzinskaya 10 123810 Moscow, Russia. E-mail: sobolev@uipe-ras.scgis.ru
Yu. S. Tyupkin
Geophysical Center, Molodezhnaya 3, 117296 Moscow, Russia. E-mail: tyupkin@wdcb.ru
See also: Sobolev, G.A., Yu.S. Tyupkin, 2000, An analysis of the
seismic energy emission process during formation of the main rupture in
laboratory experiments on destruction of rocks and before large earthquakes,
Fizika Zemly, N2, 44 - 55 (in Russian)
Microsoft INTERNET EXPLORER 5 and higher is
recommended
Contents
- Introduction
- Some grounds of the time-to-failure mode
- Analysis of the laboratory experiment data
- A retrospective analysis of strong earthquakes
- Discussion of results
- Conclusions
1. Introduction
One of the possible scenarios of the process of the main shock preparation lately discussed (Allegre et al., 1998; Ito and Matsuizaki, 1990; Keilis-Borok, 1990) is based on its analogy with the process of phase transition. An earthquake is considered as a critical phenomenon, which occurs when fracturing becomes coherently self-organized at different scales. This process develops “from below upwards” following the energy scales of self-organized fractures and is eventually concentrated in the vicinity of the hypocenter of the main shock. The laboratory experiments confirm the existence of the process of hierarchical formation of the microcracks system prior to the trunk rupture and microcracks concentration in the area of the trunk rupture at the last stage of destruction of the deformed rock specimens. (Sobolev et al., 1995).
Varnes (1989) suggested the following equation for modelling the process of energy release during the large earthquake preparation (the time-to-failure model):
(1a)
or in the integral form:
(1b)
Here k, 0<m<1, and the predicted time of the main shock tf
are unknown parameters of the failure function Q(t) (t
tf). These
parameters are defined by fitting the failure function to the time series of
accumulated Benioff strain Qexp(t)=å (Ei)1/2, which is computed from
regional seismicity in the areas that may produce a large earthquake.
Equation (1) is a special case of a more general equation, which was formulated by analogy with the phase transition theory (Sornette and Sammis; 1995):
(2)
(the log-periodic time-to-failure model).
In the next section we discuss a possible relation of the models (1) and (2) to the hypothesis of fractal structure of seismicity. The results of application of the log-periodic time-to-failure model to the analysis of the process of acceleration of seismic energy emission in the laboratory experiments on rock destruction and before strong earthquakes on Kamchatka (Sobolev and Tyupkin, 2000) and in Italy are presented in the sections 3 and 4 respectively.
2. Some grounds of the time-to-failure mode
Let us consider the following simplified scenario of the main shock
preparation. As mentioned above, the state of stress approaching the limit of
the long-term durability of crustal rocks is formed on a vast area W in
the course of preparation of a large earthquake. At first, the temporal and spatial
distribution of seismic activity in W area is more or less uniform. At
that stage, the inflow of elastic energy into this area exceeds its outflow
owing to the occurrence of weak earthquakes, creep, etc. The crustal
inhomogeneity generates a more intensive concentration of elastic energy in
some of the zones of W. As a result, the stress approaches the critical
value and the intensity of seismic energy emission (the intensity of
destruction of blocks with a low durability) increases in these zones. At a
certain moment, the configurations of the stress and strain fields define the
3D area Weff
W of the future main shock location. The critical
level of concentration of fissures is reached in this zone, and the process of
their self-organization in the vicinity of the focus of the future main shock
begins to develop on the “energetic” scales, which increase with the approach
to the moment of the main shock. This process is sustained by the
redistribution of energy accumulated in the W area into the zone Weff
(t) in accordance with the energy balance relation:
Jin(t )*Seff(t )=Jout (t )*Veff (t)
+ R(t) (3)
Here Jin(t ) is the surface density of elastic energy flow “pumped” into the Weff area, Jout is the volumetric density of the total energy flow emitted as a result of the earthquakes that occurred inside Weff, Seff is the surface of the Weff area, Veff is the volume of a set of hypocenters of the potential sources of earthquakes in the Weff area, and R(t), which is in proportion to the volume V(t) of the Weff (t) area, is the dissipation of energy inside Weff, which is not related to the earthquakes. For simplicity, we assume R(t)=e V(t).
Let us suggest three speculations.
The first is that
(4)
Here Q(t) is the failure function and a and g are the unknown parameters.
Let us introduce the characteristic linear size Leff of the Weff(t) area. Then Seff (t) ~ (Leff(t))2, Veff (t) ~ (Leff(t))d, where d£ 3 is the dimension of a set of hypocenters of the potential earthquake sources in Weff, V(t) ~ (Leff(t))3. The equation (4) may be re-written as follows:
(5)
It is natural to suppose that Jin(tf)¹ 0. Then Jin(t)= J0+ J1 (tf-t)Tc-1+…when time t ® tf . (Here Tc is the characteristic time of the preparation process of the main shock.).
The second speculation is that in the course of preparation of a strong earthquake, the process of destruction of the medium collapses towards the place of the future main shock (see, for example, Zavyalov and Nikitin, 1999. )
Finally, the third speculation is that the Leff(t) function is analytical at the critical point tf. Then
when
t ® tf (6)
As a result, the integration of equation (5) produces, in the case if t ~
tf
(7)
Here and below O(z) indicates the members tending to zero more quickly than z at z® 0 and
.
If 2<d<2+a-1,
then equation (7) coincides with the Varnes equation (1b) with m=1+(2-d)a.
The speculation (6) can be reformulated as the following:
for t® tf (6a)
where
and l is arbitrary.
The assumption that Leff(t) is an analytical function in the vicinity of the critical point tf may not be correct since the scenario is developed on a fractal manifold. However, the reasons discussed above about the scaling invariance of the process at the stage, when the state of stress is close to the limit of long-term durability of crustal rocks, allow, by analogy with the theory of phase transition, to replace the assumption (6a) by the assumption that
Leff(t)=l0-1
Leff(Sl0(t))
for t® tf (6b)
where Sl0 (t) is the discrete scaling transformation Sl0 (t): (t-tf)® l0(t-tf) (Sornette 1998). In this case, a general expression of the singular part of a function F(t) in the vicinity of the critical point tf will be like
,
where z=|t-tf|, H(z) is the periodic function: H(z+1)=H(z), and b is defined by the linear part of the Sb (t)) transformation in the vicinity of the critical point tf : Sb (t))»tf +b(t- tf ) and d is an arbitrary constant (see, for example, Derrida et al. 1983). Therefore, the replacement of the Leff(t) function analyticity assumption at the critical point tf by the assumption that this function is invariant under the renormalization group transformation Sb , produces the following generalization of (6) in case of t ~ tf :
(8)
where h(1+z)=h(z). By referring to the results of application of the model to the experimental data, which is discussed below, we consider co<<1. We also use the principle that the minimal number of formal parameters, which have no clear physical interpretation, must be used in the model and d=1 is selected. Introducing (8) into (5) and approximating the periodic function h by a cosine, after integration of equation (5) we obtain:
(9)
where c*~ c0<<1, and parameters p
and f are expressed
in terms of parameters b and
Tc. The expression (9) changes into (7), if c0=0.
The parameters of the failure function Q(t) are defined by fitting Q(t)
to the time series of accumulated Benioff strain Qexp(t)=å (Ei)1/2. Here the
sum involves the earthquakes that occurred in the vicinity of the hypocenter of
the main shock before time t; Ei is the seismic energy of
these earthquakes. It is convenient to present Qtotal=Q(tf)
as Q(t1)+A1, if we have the information about
seismicity in the vicinity of the main shock till time t1< tf.
In this case, the Q(t1) is known, and A1 is
generally determined by the energy of the main shock.
To compare the model (9) with the experimental data, we used the
following expression of function Q(t):
(10)
and the coefficients A1, A2, m, tf, c, p, f are determined by minimization of quadratic form
(11)
(by the least-square method)
At the end of this part, let us note that the independent assessments of the
power m and the fractal dimension d give the estimation of the
parameter a, which was
introduced in (4), and the relation between the seismic energy and the total
energy released as a result of an earthquake is nonlinear in the general case.
3. Analysis of the laboratory experiment data
The power m exponent in (10) is the function of d, which is interpreted as the fractal dimensionality of manifold of potential hypocenters. This means that quantity m can vary according to the degree of fracturing of the medium. To test this hypothesis we applied model (10) to the laboratory research data. In a number of papers the authors indicate a similarity of statistical properties of the acoustic emission process in laboratory rock specimens subjected to high pressure and seismicity (see, for example, Sobolev, 1993). In particular, it was established that the distribution of the acoustic emission impulses by energy corresponds to the Gutenberg-Richter relation of frequency-magnitude statistics of earthquakes, it was also found that the spatial distribution of the acoustic emission sources has the fractal structure typical of seismicity as well, and so on. In this part of the paper we present the results of application of equation (10) for modelling the process of energy release as a result of acoustic emission, which was noted in the course of laboratory experiments AE36 and AE 42 (Lockner, Byerlee, 1980; Lockner et al., 1992, Ponomarev et al., 1997). By virtue of the conditions of the experiment, the different stages of destruction were studied on specimen AE42 and on specimen AE36. In fact, for specimen AE42 the experiment was ended after passing through the peak load at the level of 80% of the maximal stress. The acoustic emission data of the experiment reflect mainly the initial phase of formation of the main rupture. In experiment AE36, the first burst of acoustic emission occurred in the vicinity of 25000 s. In the 31500 - 32500 s. interval, with the load of about 55% of the maximal, the rupture reached the lower base, while the rupture plane reached the middle of the bearing piston of the specimen. This process initiated another increase of stress. The load grew up to 37500 s. and caused the formation of a fresh plumage of macroruptures in the lower part of the specimen followed by an intensive acoustic emission that started after 43368 s. at the final stage of the experiment. Moreover, the growth of the rupture in this specimen was rather irregular, with acceleration and slowing down as a consequence of destruction of the minor inhomogeneities of the medium. The plots of acoustic energy emission are given in Figure 1(Figure 1a - experiment AE36 and Figure 1b - experiment AE42). In this case, the acoustic energy was summed up in the intervals of 110 s and 10 s, respectively.
For the AE36 specimen, an analysis was carried out of the acceleration intervals in the acoustic emission process at the moment of the first (23075 - 25275 s) and of the last (40205 - 44255 s) intensive burst of acoustic energy emission. For the AE42 specimen, the 7664 - 15364 s. time interval was studied, which covered the initial stage of formation of the trunk rupture. (The peculiar feature of these experiments was that the growth of acoustic emission rate caused retardation of the axial deformation in order to avoid a rapid destruction of the specimen. Therefore, in contrast to real seismicity, the tf value in these experiments cannot be interpreted as the time of the trunk rupture).
A large amount of data obtained as a result of the experiments AE36 and AE42 allows to estimate the fractal dimensionality of the manifold of acoustic emission sources, which were recorded at the intervals, when the process of acoustic energy emission is accelerated, and to compare these estimations with the corresponding values of parameter m in model (10), thus providing a possibility to obtain the value of the a parameter introduced into (4).
Figure 2 shows the plots of failure function Q(t)
for the values of parameters having the best correlation with the experimental
values Qexp(t)=å
(Ei)1/2 of the acoustic emission data in the
experiment AE36 (the time intervals are 23075 - 25275 s.,
Figure 2a ;
and 40205 - 44255 s., Figure 2b ) and in the
experiment AE42 (time interval 7664 - 15364 s., Figure 2c ).
Here Ei is the energy of the i-impulse of the sequence
of events that occurred by the moment of time t in the 
3 mm. layer near
the plane of the rupture. The experimental values Qexp(t)
were calculated with the time interval of 50 s. The obtained parameters of the
model are given in
Table 1. The inset shows the plot
of the
4. A retrospective analysis of strong earthquakes
We have carried out a retrospective analysis of strong earthquakes that
occurred on Kamchatka and in Italy. For Kamchatka, the data of the Kamchatka
instrumental catalog compiled by the Kamchatka Expedition for Experiments and
Methodology of the Institute of Volcanology, Far East Department of RAS (1962-1995),
supplemented by the data of the operative catalog (1996-1998) was used. The
energetic class K is presented in this catalog instead of the magnitude M as
the energetic characteristic of an earthquake. The catalog presents,
practically without gaps, a list of earthquakes of magnitude M
3 (energy class K9) that occurred
within the area of the eastern coast of the central and southern Kamchatka (Fedotov, 1987). The events with M7.0 (K14.5) and with the depth of the
hypocenter <90 km were analyzed for the period of 1965-1998.
Table 2 shows parameters of these earthquakes. The
magnitudes published by “Earthquake in the USSR” issues are pointed out in this
Table. One can see that the estimation of energetic class and magnitude is a
little bit differs for some earthquakes. Therefore the energetic class is used
as an energetic characteristics when we present the results of analysis of
Kamchatka earthquakes (see Figure 3 and
Table 4). An asterisk indicates earthquakes with
epicenters located in the area apparently not covered by reliable data on weak
events in the catalog; therefore, these earthquakes were not analyzed.
For Italy, the earthquakes with magnitude M>5.0, which are listed in the Table 3, were analyzed. These earthquakes were chosen because their epicenters are located in the regions where weak events are recorded with sufficient reliability. We used the seismic catalog compiled by the Istituto Nazionale di Geofisica, (ING) for 1975-1998 period (Barba et al., 1995). The completeness of the ING catalog is defined by the magnitude Mmin=3.5 for 1975-1985 (Wyss et al., 1997) and by the magnitude Mmin=2.3 for 1986-1998. Both catalogs were cleared of aftershocks according to the algorithm of Molchan and Dmitrieva (1991).
When using model (10) in the analysis of seismicity, we should bear in mind the following parameters:
- Parameters of analyzed seismicity.
·
The characteristic linear dimension of the analyzed
region Lc. We have chosen Lc =1.5
Reff, where
Reff is the characteristic linear dimension of the source of
predicted event determined by its energy characteristics. In the case of
Kamchatka, the empirical relation was used between Reff and energetic
class K, which has the form of log(Reff)=0.244K-2.266
(Riznechenko, 1976). For Italy the relation log(Reff)=0.44M-1.289
was applied (Papadopoulos and Voidomatis,
1987)
·
The depth interval (Zmin, Zmax)
of events preceding the main shock, which are applied for the analysis. The
choice of Zmin, depends, generally speaking, on the depth of
the focus of predicted event. It was shown (Sobolev and
Tyupkin 1997, 1999) that on Kamchatka the process of seismic quiescence
before large earthquakes with the hypocenter depths 
30 km, which later changes to
foreshock activation, is best observed by the analysis of weak seismicity at
more than 20 km depths of hypocenters of events. For Kamchatka we used Zmin=0
in the analysis of the surface event of 1997 and of the earthquake of 1971 with
the depth of the hypocenter at 25 km; Zmin=20 km is applied
for other events. As for the Zmax value, its
choice was determined by the circumstance that for Kamchatka the intensity of
weak seismicity flow quickly attenuates at the depth more than 70-80 km. We,
therefore, used Zmax=100 km, and for the case of the surface
earthquake of 1997 it was verified that the result does not depend on the
choice of Zmax value within the 50-100 km interval.As shown
in
Table 3, the depth of the hypocenters of the
analyzed earthquakes of Italy is not more than 20 km. Therefore, in this case,
we used Zmin=0 km and Zmax=100 km.
· The interval of magnitude of events preceding the main shock that is used for analysis. The level of completeness of the catalog determines the minimal magnitude Mmin. As the maximal magnitude, we did not impose on its value any restrictions.
·
The characteristic time interval T0 of
the process of self-organization of fissures in the vicinity of the focus of the
future main shock, which depends on its magnitude. It was assumed for Kamchatka
that T0 is about 4 years for the earthquakes with 7
M<7.4
and 5 - 8 years for the earthquakes with 7.4M; moreover, the T0
variations within two years of extension of these intervals had no significant
effect on the result. We chose T0 ~ 2 - 3 years for the
analyzed earthquakes of Italy.
· Maximal duration of temporal interval of prediction DT. Owing to the limited number of available strong earthquakes for analysis, we did not study the possible estimation of DT value.
- The m degree exponent which is the parameter characterizing the “ fracturing” of the medium in the vicinity of the future earthquake focus.
- The “trend” model parameters:
· Parameter A1. It determines the predicted magnitude Mmod (energetic class Kmod) of the main shock. For Kamchatka the energetic class is used as the energetic characteristic of an earthquake and Kmod=2 log(A1); for Italy Mmod=1.39 log(A1) - 3.67. (We use the following relation between the seismic energy E of earthquake (in Joule) and its magnitude M for Italy: log(E)=1.44 M+5.29 (Papadopoulos and Voidomatis, 1987)).
· The predicted time tf of the main shock.
· Parameter A2.
- Parameters of the log-periodic
fluctuation: amplitude c, phase f,
period of log-periodicity p
As discussed above, the m degree exponent characterizes the “fracturing” of the medium in the area of the focus of the future main shock. If we suppose that the statistic properties of “fracturing” of the medium of the region were formed during a long time and can not be changed in the course of several tens of years, the m parameter can be determined by the “training” events. We used the earthquake of March 2, 1992 and the Umbria-Marche earthquake of September 1997 as the “training” examples for Kamchatka and for Italy, because the “acceleration” process was clearly manifested before these events. In order to determine the unknown coefficients of the model (10), the least square method was applied; i.e. the expression (11) was minimized by all seven parameters. In the case of “training” events, the main shock also participates in the determination of the coefficients of the model. In both instances the value of the m degree exponent was obtained = 0.12.
For the remaining earthquakes, the parameters A1, A2, tf, c, p, f, were determined from minimization of the expression (11) at the fixed value of parameter m=0.12 and under condition that the main shock does not participate in (11). We made the same estimation of parameters of the model (10) for the “training” earthquakes as well. We guessed that the model gives a wrong prediction if the difference of the predicted time tf and the real time tr of main shock is more than one year.
The results of “prediction” of the model (10) for Kamchatka earthquakes are shown in Table 4 and in the Figure 3, and for the earthquakes in Italy in Table 5 and in Figure 4. The model gives wrong prediction (tf - tr > 1 year) for the Kamchatka earthquake of 28/02/1973 and for the Italian earthquake of 23/11/1980
5. Discussion of results
As prognostic, the model has obvious drawbacks. First, it can in no way determine the place of the supposed earthquake and, consequently, the area that is to be analyzed with its help. Second, the model does not distinguish earthquake swarms from foreshock activation. Figure 5 shows the results of the model’s “prediction” for the point with coordinates Lat=41.5oN, Long=14.6oE, in the vicinity of which a rather intensive swarm was observed in the middle of September 1996. But the model can help with the estimation of the magnitude of the expected event and in obtaining a more precise moment of the predicted main shock, if the model is considered as a supplement to the prognostic approach, which determines the potentially hazardous area and the time interval of the alarm.
One of the problems in the application of this model is that it implies a certain degree of isolation of the analyzed area from the impact of other strong earthquakes during the T0 time interval. In the case of Kamchatka earthquakes that occurred on July 8 and November 13, 1993 at a distance of 96 km between them, this condition is obviously not fulfilled, and the space-time proximity of these events could have caused the peculiar features in the development of the seismic process in the vicinity of their epicenters.
The last comment is associated with the analysis of the Kamchatka earthquake that took place on December 5, 1997. The regional Kamchatka seismological network, the Obninsk Seismological Station and the World Seismological Network reported considerably different coordinates of the epicenter of this earthquake. In the present study, the determination of the coefficients of the model (10) was made by using data on earthquakes that occurred in the radius of 37 km from the center located at the point of maximum foreshock activation observed before the main shock. The map of spatial distribution of foreshock activation is presented in the Figure 6 (Sobolev and Tyupkin, 1999). An asterisk shows the location of the instrumental epicenter of the earthquake with coordinates Lat=54.95oN, Long=163.29oE reported by the regional Kamchatka Seismological Network. It is apparent that it is located near the center of foreshock activation. Since for this earthquake the acceleration of the process of seismic energy release before the main shock is demonstrated with sufficient clarity, we made an additional estimation of the parameters of the model by minimizing the quadratic form (11) in all seven parameters A1, A2, tf, c, p, f, m. Accordingly, the value of power m index was found to be equal to 0.12, the same as in the case of the earthquake of March 2, 1992.
6. Conclusions
The analysis of the process of acceleration in the release of energy in the course of destruction of rock specimens and before strong earthquakes on Kamchatka and in Italy shows that:
- There are grounds to believe that the value of power m index in the model (10) is associated with the dimensionality of the potential earthquake hypocenters manifold in the vicinity of the epicenter of the future strong earthquake.
- Model (10) approximates well the acceleration of the process of acoustic energy release in the course of development of the trunk-line rupture in experiments AE36 and AE42 on destruction of rock specimens. During these experiments, at different stages of development of destruction, the parameter m of model (10) changes from 0.09 to 0.26. An independent estimation of correlation dimension d of the manifold of acoustic emission sources at the acceleration stages of the process of acoustic energy release correlates with these changes in parameter m. The comparison of m and d values allows to estimate the value of parameter a introduced into (4), which defines the law of changing of the seismic and total energy release ratio of an earthquake as the function of its magnitude.
- The acceleration of seismic energy release is not a universal process that is observed before all strong earthquakes. For example, this process is detected in the earthquakes that occurred on the eastern coast of the Central and Southern Kamchatka mainly at the depth of the hypocenter <50 km. The application of model (10) to estimate the magnitude and moment of the future strong earthquake is possible in the case, when a tendency to acceleration is observed in the release of seismic energy in the vicinity of its epicenter.
- There are grounds to believe that the parameter m in the model (10), in the first approximation, can be considered universal for those areas of the region, where the acceleration of seismic energy release is observed before strong events (see also Bufe and Varnes, 1993). According to our estimation, m=0.12 for the eastern coast of the Central and Southern Kamchatka and for the Central Italy.
- In all the cases analyzed herewith the log-periodic component exists, though the amplitude c of this component is always <<1.
- As prognostic, the model (10) has several drawbacks. It does not determine the place of the supposed earthquake and, therefore, the area, which should be analyzed with its help. Moreover, it does not distinguish between the earthquake swarm and the foreshock activation. But if we consider this model as supplementary to prognostic approach, which determines the potentially hazardous areas and the time interval of alarm, it can help to estimate the magnitude of the expected event and to more precisely determine the prognostic moment of the main shock.