Tectonophysics,
26 (1975)
0 Elsevier Scientific
91-101 Publishing Company, Amsterdam - Printed in The Netherlands
A STOCHASTIC MODEL OF EARTHQUAKE OCCURRENCE AND THE ACCOMPANYING HORIZONTAL LAND DEFORMATION
YUKIO HAGIWARA Earthquake
Research
Institute,
University
of Tokyo,
Tokyo
(Japan)
(Submitted
June 19, 1974; revised version accepted October 29, 1974)
ABSTRACT Hagiwara, Y., 1975. A stochastic model of earthquake occurrence horizontal land deformation. Tectonophysics, 26: 91-101.
and the accompanying
A large-scale earthquake is believed to be associated with a release of strain energy accumulated in the crust, probably by the motion of upper-mantle lithosphere. Such an earthquake mechanism is well simulated by a belt-conveyer model proposed by Utsu (1972). The probability of earthquake occurrence can be estimated on the assumption that the motion of a slider on the belt-conveyer is mathematically formulated as a Markov process. In the probabilistic expressions, the results of Mogi’s (1962) rock-fracture experiments are applied to the hazard-rate function of earthquake occurrence. The hazard-rate function has two coefficients, A and B, to be determined by the experiments. It is concluded that, when B is small, a number of small-scale earthquakes occur in the early time after the accumulation of crustal strain energy starts, but that the accumulated strain energy changes catastrophically into a single large-scale earthquake, when B is large.
INTRODUCTION
It is the ultimate goal of earthquake-prediction research to foretell the time and place of occurrence and the magnitude of predicted destructive earthquakes. Place of occurrence and magnitude can .be estimated to some extent by investigating local features of earthquake occurrences from records and experiences in the past. Observations of microseismic activities, measurements of changes in seismic-wave speeds, repeatedly made geodetic surveys, detections of crustal tilt, strain, specific resistivity, etc. are thought to be useful in tackling the occurrence-time problem. The occurrence is, however, governed by so many physically uncertain elements that it is extremely difficult to forecast the time when a large-scale earthquake will occur. Despite the difficulty in predicting occurrence-time, a quantification of earthquake warning is urgently required for security against earthquake damage. Rikitake (1969) presented an attempt to apply probabilities of occur-
92
Fig. 1. Utsu’s model. The belt-conveyer corresponds to a lithospheric plate moving towards a continent (the wall in this figure). An elastic motion of the spring becomes an earthquake.
rence-time and magnitude to such a problem. It may be certain that such probability considerations are one of the basic approaches to a long-range forecast of earthquake occurrence, on the basis of which personnel distributions and investments in equipment can be properly controlled to minimize disastrous damage. Furthermore, Rikitake (1974) proposed a method by which the probability of occurrence-time can be estimated on the basis of the statistics of ultimate strain of the crust as obtained from geodetic surveys repeatedly made over earthquake areas. Hagiwara (197413) later re-examined the same data by means of a Weibull distribution analysis. According to the current theory of plate tectonics, the crust of northeastern Japan is strained by a westward drifting of the Pacific Plate of the lithosphere, and similarly the crust of southwestern Japan is strained by the motion of the Philippine Sea Plate. An earthquake occurring off or along the Pacific coast of the Japan Islands can be explained as an elastic rebound of the strained crust against friction acting on the crust-plate boundary surface. The strain energy accumulated in the crust is released partially or totally by the earthquake. Utsu (1972) put forward a simple model simulating such a type of earthquake mechanism. The Utsu model is composed of a belt-conveyer, a weight and a spring as shown in Fig. 1. The lithospheric plate corresponds to the belt of the model, and Japan’s crust to the spring connecting the weight with a stable wall (continent). As the bottom surface of the weight is slightly rough, it is conveyed on the moving belt, but suddenly it slides back due to elastic reaction of the spring. When the weight is pushed back, strain energy accumulated in the spring turns into an earthquake and the associated spring extension corresponds to the horizontal land deformation. It is thought that the motion of the weight in the Utsu model may be mathematically formulated as one of the problems which are amenable to analysis of the Markov stochastic process. On the analogy of the Utsu model, the stochastic behavior of horizontal land deformations can also be described
93
in the form of the Markov transition probability. The present paper contains the probability expressions for the Utsu model using the hazard-rate function, the number of fracture occurrences per unit time, which Magi (1962) obtained from the results of rock-breaking tests. TRANSITION
MATRIX
EXPRESSING
HORIZONTAL
LAND
DEFORMATION
In order to formulate the stochastic behavior of horizontal land deformations associated with the motion of the lithospheric plate and earthquakes, the emphasis is on a probabilistic treatment of the Utsu model. Suppose an Utsu model having a number of weights on the belt-conveyer. Fig. 2 shows an example of the model having three weights A, I3 and C, which are connected with the same wall and move independently on the belt-conveyer. Assuming that strain .c takes on discrete values designated by n = 1, 2, . ..N. we define: E = (n-l)Ac
(1)
where Ac is a strain interval. As seen in Fig. 2, A and C are located at n = 2, but B at II * 1. Strain increasing at a rate u during a short time interval At, we have: Ae = uAt
(2)
Let p,(t) be the probability of the event that a weight is located time t. The corresponding probability vector is given by: PI(t)
P(t) =
Pz(t) :
1
where elements
of the probability
at n at
(3)
vector,
in general, have to satisfy a relation:
N
c p,(t)
=
1
n=l
In Fig. 2, we see that pl(t) = l/3, p2(t) = 2/3 and p3(t) = p4(t) = 0. Assuming that no earthquake occurs in a time interval between t and t + At, the weights shift their locations from n to n + 1 on the moving belt. In the case of the example shown in Fig. 2, we obtain p2(t + At) = l/3, p3(t + At) = 2/3 and pl(t + At) = p4(t + At) = 0. If we express the relation between the probability vectors at t and t + At as: P(t + At) = W(t)
(5)
Fig. 2. Utsu’s model with three weights (sliders).
the transition
matrix R has a form:
Furthermore, we suppose another case that an earthquake accompanied with a horizontal land deformation occurs in a time interval between t and t + At. In this case, A and/or C may slide back from n = 2 to 1, although B does not move any more. When a weight shifts its location, in general, from 171to n, the condition m 2 n should be satisfied. The ~o~esponding t~nsi~on matrix is given by:
(7)
where the eIements of Q are non-negative. course to be satisfied:
E n=l
Qmn =1
The following
condition
has of
(8)
95
The relation between the probability vectors at t and t + At has a form similar to eq. 5: P(t + At) = QRP(t)
(9)
In the right-hand side of the above relation, the transition matrix from t to t + At is not Q but the product of Q by R. This indicates that the crust is continuously strained by the moving plate even at the time when an earthquake occurs. STOCHASTIC
PROCESS
OF HORIZONTAL
LAND
DEFORMATION
The probability of occurrence of an earthquake accompanied by horizontal land deformation in a time interval between t and t + At being denoted by AAt, the probability of the event that such an earthquake does not occur in the same time interval is 1 - hAt. We call X the “hazard rate” in terms of probabilistic quality-control techniques. The stochastic processes eqs. 5 and 9 are assigned to the events with the probabilities hAt and 1 - AAt, respectively. These two events are mutually exclusive, so that, combining eq. 5 with eq. 9, we can finally obtain the equation governing the stochastic process of the motion of the weights of the Utsu model in the form: P(t + At) = [(l - hAt)I + hAtQ]RP(t)
(10)
where I is a unit matrix. By using elements of the matrix, this formula can be written as:
p,(t
+ At) =
(1- A,-&lp,-,(t)
+XNqN,$N(t)
PN(t + At) = [I-
+ [1 - b&(1
+ At
1
c L-a,znPm-l(t) 1??l=?l
(11)
for2Sn
&v-&(1 -
-
qNN)bN-l(t)
qNN)bN(t)
where X is replaced by X, because h depends on n. Starting with an initial condition:
11 1 0
P(0) =
.
(12)
96
and using eqs. 11, we can numerically determine all the values of p,, (t). In order to solve eqs. 11 practically, however, we have to give the hazard function and the elements of the matrix Q. The hazard function is derived from Mogi’s (1962) empirical formula based on results of rock-breaking tests. In has a form: X = A exp(Bo)
(13)
where A and B are constants to be empirically determined and u is the external stress. On the analogy of the fracture occurrence in a rock specimen, the above relation may also possibly hold good for earthquake phenomena. Rikitake (1974) calculated the ultimate strain of the earth’s crust from the geodetic data of repetition surveys made before and after earthquake occurrences. Based on Rikitake’s table, Hagiwara (1974a) concluded that an exponentially increasing hazard model fitted well with the statistical distribution of ultimate strain, and that the value of B determined from the ultimatestrain statistics was almost equal to Mogi’s value determined from laboratory experiments. This may imply that B is independent of the size effect. If A is obtained from the geodetic measurements, the same hazard function as eq. 13 can be applied to the crustal fractures without any size-effect considerations. Assuming that Hooke’s law is satisfied as a stress--strain relation, we have: O=EE
(14)
where E is Young’s modulus of the earth’s crust. Substituting 14 and taking eqs. 2 and 13 into account, we finally obtain: X, = A exp[BEAe(n Moreover,
eq. 1 into eq.
- l)]
we assume that, if we know B and N, A is determined
AAt = exp[-BEAe(N
- l)]
(15) from (16)
On the other hand, the elements of the matrix Q can be somewhat arbitrarily defined as long as the condition (eq. 8) is satisfied. As a matter of fact, the formula of Q is unknown. It may not be theoretically derived from the Gutenberg-Richter relation together with some assumptions of strain-magnitude relation, because the motion of the weight of the Utsu model involves theoretically uncertain probabilistic behaviors due to friction between the weight and the belt. Instead of sophisticated procedures for theoretically determining the probability, we adopt a simple and convenient distribution which well describes the problem. If a binomial distribution function of a variable ~(0 < /J < 1) is adopted, an element of the matrix Q becomes: m-l
(I-
4 mn = I n-l
p)n-lpm-n
(17)
I
We see that eq. 17 satisfies eq. 8. It is convenient for practical to use, instead of eq. 17, the following recursive formulas:
computations
97
qn&n+1=
(m - n)(l - PJq Inn
nP
Provided that the probability vector has a stationary solution, it is obtained by puttingp,(t + At) = p,(t) in eq. 11 and solving the simultaneous equations of N variables. In this case, we adopt the relation (eq. 4), instead of the first or last equation of formula 11. COMPUTATION
RESULTS
In order to obtain the stochastic process of horizontal land deformation, At, B, N and /L have to be given. If u is known, we can determine Ae from At through eq. 2 and A from At, B and N through eq. 16. In the present study, we take At = 2 years and u = 5 - lo-‘/year, so that Ae = 10H6. Furthermore, we take N = 50. B and p, which characterize the process, are taken as parameters. Mogi (1962) obtained A = 4.0 - 10-18/min. and B = 0.37 cm2/kg in laboratory experiments of granite fracture. A depends on a size of rock sample, so that Mogi’s value of A is not applicable to the crustal problems. According to Hagiwara (1974a), however, B is independent of a size effect, so that Mogi’s value of B may be used for our present purpose. Assuming Young’s modulus asE=5-105k g / cm2 and referring to Mogi’s value of B, we take B = 0.2, 0.4 and 0.8 cm2/kg in the present computations. Figs. 3A, 3B, 3C and 3D show the stochastic processes starting with the initial condition (eq. 12) for the cases of B = 0.2 cm2/kg and p = 0.1, B = 0.8 cm2/kg and ~1= 0.1, B = 0.2 cm2/kg and p = 0.9, and B = 0.8 cm2/kg and 1-1= 0.9, respectively. It may be seen in these figures that, if B is small, a number of small-scale land deformations occur in the early time after the crustal strain energy accumulation starts. In this case, the land deformations may be associated with many small-scale earthquakes like an earthquake swarm. On the other hand, a land deformation with large B appears catastrophically when the accumulated strain energy changes into a single large-scale earthquake. p is also related to the scale of the earthquake. Small-scale earthquakes occur when p is small, but large-scale ones occur when /.Lis large. A land deformation associated with a large-scale eartliquake occurring off the Pacific coast of Japan can be simulated by the model in Fig. 3D (B = 0.8 cm2/kg and /l = 0.9). One of the most interesting matters to be studied is the behavior of horizontal land deformations associated with a catastrophic large-scale earthquake. We assume an event that the weight of the Utsu model shifts its position from m to n in the time interval t - t + At. If n is taken as n = m - K, K is an integer indicating the length of the slide, that is, the scale of land deformation.
98
B = 02cm2kg-' P = 01
4 jI
I
0 20-
co 40 z ? z 60 s ; 80
100
m-
5
4
3
2
I
x 12
STRAIN
A
.E 0
6 = 08cm2k$ y = 01
I
oi-----5 0
4
3
2 STRAIN
I
0 x105
99
B = 02 ,m2 kg-' )I= 09
m
5
4
2
3
C
1
STRAIN
----E 0 x Id'
B = 08cm'krj' p=
m5 Q
09
_E 4
3
2 STRAIN
1
0 xt$
Fig. 3. Probability p,(t) in the cases of: (3A) I? = 0.2 cm2/k#, p = 0.1; (3B) B = 0.8 cm2/kg, /J = 0.1; (3C) B = 0.2 cm2/kg, I_(= 0.9; and (3D) B = 0.8 cm /kg, p = 0.9. The initial values (t = 0) are taken as 1 at E = 0.
TIME
Fig. 4. Probability (K = 20).
IN
YEARS
GK(t) assuming
that /J = 0.9 and crustal
strain is greater
The probability of the event that the land deformation than K is given by: N-K @K(t)
= At
c
than 2 * lo+
has a scale greater
N c
&,&n&(t)
m=K+lm=n+K
= At
5
h&,(t)[l
-B,Jm
-K,
K)/B(m -K,
K)]
(19)
m=K+l
where B and B, are the beta and incomplete beta functions. Fig. 4 shows the probability @K(t) of the cases of B = 0.2; 0.4 and 0.8 cm2/kg, assuming /J = 0.9 and K = 20 (i.e., the released strain is larger than 2 * 10m5). It is seen in the figure that, when B is small, deviations from the mean occurrence time are relatively large, When B is large, the occurrence time becomes almost cyclic. That is to say, it can be relatively easy to foretell the earthquake-occurrence time in an area of large B on the basis of observation data of crustal strain. SUMMARY
The earth’s crust is horizontally deformed as strained by the convective motion of the lithospheric plates. The behavior of the crustal motion is well simulated by the LJtsu (1972) model. The hazard rate, which Mogi (1962) obtained from laboratory experiments of rock fracture, is used for the formulation of the Utsu model as one of the problems available for the Markov stochastic process. The elements of the transition matrix are assumed to be binomiaI,distribution functions having a parameter p, which is related to the scale of land deformation. A small-scale land deformation occurs when p is small, but a large-scale one corresponds to a large value of p.
101
The Mogi hazard rate has two coefficients, A and B, to be empirically determined. The value of B is independent from a size effect (Hagiwara, 1974a), so that it is applicable to statistical problems of ultimate strain of the crust. We conclude from our computation results that, when B is small, a number of small-scale land deformations occur in the early time after the crustal strain-energy accumulation starts, but when B is large, the accumulated strain energy changes catastrophically into a single large-scale earthquake. This conclusion can also be qualitatively anticipated. As a result, combinations of p with B make a number of probabilistic patterns of horizontal land deformation. If 1-1and B are statistically determined from the past records of land deformations and results of rock-breaking tests, we can estimate to some extent the occurrence time and magnitude of a large-scale earthquake expected to occur, which is accompanied by a land deformation. Generally speaking, a probabilistic interpretation is necessary only because of our ignorance of a physical law governing the deterministic evolution of nature. The theoretical development in the present study rests upon a certain statistically assumed mechanism for earthquakes. If a physical model of fracture mechanism is established in both the crust and a rock piece, results of rock-breaking tests can be applied to earthquake phenomena without physical uncertainties in the analogy between nature and laboratory. Fracturing tests from small to giant-sized rock pieces may perhaps be the easiest way to solve such a problem as to what extent laboratory results can be transferred to the earth. ACKNOWLEDGEMENT
The author wishes to express his cordial thanks to Professor T. Rikitake. It was at his suggestion that the author set about the present study.
REFERENCES Hagiwara, Y., 1974a. Probability of earthquake occurrence estimated from results of rock fracture experiments. Tectonophysics, 23: 99-103. Hagiwara, Y., 1974b. Probability of earthquake occurrence as obtained from a Weibull distribution analysis of crustal strain. Tectonophysics, 23: 313-318. Mogi, K., 1962. Study of elastic shocks caused by the fracture of heterogeneous materials and its relations to earthquake phenomena. Bull. Earthquake Res. Inst., Univ. Tokyo, 40: 125-173. Rikitake, T., 1969. An approach to prediction of magnitude and occurrence time of earthquakes. Tectonophysics, 8: 81-95. Rikitake, T., 1974. Probability of earthquake occurrence as estimated from crustal strain. Tectonophysics, 23: 299-312. Utsu, T., 1972. Aftershocks and earthquake statistics, IV. J. Fat. Hokkaido Univ., 4: l-42.