Geophys. J. R. astr. Soc.

(1987) 88, 723-731

Random stress and earthquake statistics: time

dependence

Y. Y. Kagan and L. Knopoff lnstituteofGeophysicsandP!anetary

Physics, University of California, Los Angeles, California 90024, USA

Accepted 1986 August 22. Received 1986 August 22; in original form April 24.

Summary. The interearthquake time distribution is analysed on the

assumption that those stresses that are not observed directly, change in a way
that is describable by a random walk, i.e. as a Brownian motion. In this case,
the time intervals between earthquake pairs has a power.law distribution with
exponent -3/2. If tectonic stress loading is added to the Brownian motion,
the interearthquake time distribution changes from a power-law to an almost
Gaussian renewal distribution. The actual distribution depends on the ratio of
the size of the random component to that of the tectonic component. We
find that after about two days, tectonic stresses influence the temporal
distribution of aftershocks, for main shocks with ML = 1.5. In the random
regime the Omori law holds, while in the tectonic regime, an exponential
distribution holds. The time of transition between random and tectonic
effects increases with the size of the main shock.

Key words: random stress, earthquake statistics, Omori law

Introduction

In earlier papers (Kagan & Knopoff 1981; Kagan 1982) we have described a kinematic
stochastic model of an earthquake process, which generates most, if not all, of the statistical
features of earthquake occurrence identified thus far. Since the model is based on a small
number of arbitrary assumptions, it follows that the model is phenomenological, even
though the number of such assumptions is small. For the purposes of this paper, only one
of these assumptions is important and may be summarized as follows: the probability of
occurrence of one infinitesimal earthquake in the wake of another such infinitesimal earth-
quake is a power-law function of time with an exponent for shallow seismicity that is close
to -3/2.
The major item of seismological evidence which supports this assumption is the well-
known Omori law of aftershock occurrence. We have shown that the above assumption leads
to the Omori law directly (Kagan & Knopoff 1980, l 986a); this law is generally valid, not
24
724 Y. Y. Kagan and L. Knopoff
only for the time distribution of aftershocks, but also for the distribution of foreshocks;
indeed it explains almost all of the temporal clustering of shallow earthquakes that is reliably
known at present.
Because of its phenomenological nature, the model cannot be extended significantly
beyond what has already been accomplished. In particular, we have no convenient way to
evaluate the influences of the steady stress increase due to plate tectonics or other long-term
effects on the statistical output from the model; while the present model adequately
reproduces short-term statistical phenomena, it cannot reproduce long-term behaviour. In
this paper we rectify this defect by introducing a physical basis for the original phenomeno-
logical model. We show that a rational modelling of the fluctuations in the stresses near a
crack due to all the other cracks in the neighbourhood will allow us to duplicate the
phenomenology of the original model. This interpretation will permit us to add tectonic
stresses to the model as well. An inversion of the procedure will then allow us to identify the
statistical consequences of the superposition of both the steady and the random processes.
We show that the exponent -3/2 in the phenomenological model arises as a natural
outcome of the conventional assumption that an earthquake occurs when the stress in the
vicinity of a fault edge exceeds some critical value. In the absence of any information about
tectonic stresses, we can assume that the stress is controlled by random factors. Our aim is to
show that the value of the exponent is not an adjustable parameter, but is the consequence
of a well-defined physically based assumption. Once we have established the physical basis,
we can proceed to the discussion of the other factors of earthquake occurrence, such as the
size distribution and spatial statistics which have been discussed phenomenologically in the
above papers. We consider the problems of the relationship between the physical and
stochastic modelling of the spatial earthquake parameters in a companion paper (Kagan &
Knopoff 1986b).

2 Interearthquake time statistics

As we have indicated, in our model (Kagan & Knopoff 1981) the time intervals between
infinitesimal events that constitute our earthquake source-time function, are distributed
according to a power-law distribution with exponent -3/2. Because of the infinity in the
distribution at t = 0, we arbitrarily impose a 'dead-time' t 0 before which no new earthquakes
can occur. A physical argument also suggests that a short characteristic time be introduced.
In our model the dead-time arises because we bypass consideration of the dynamical nature
of fracturing. If we consider a fracture to be an accelerating quasistatic process, then the
seismic event itself is the culmination of the quasistatic growth phase at the moment when
the rate of fracture growth approaches seismic wave velocities. Because we bypass the
problems of dynamical rupture theory, the radiation from catastrophic fractures is taken
to be instantaneous in the model; we introduce a very small time constant to describe the
finiteness of the rupture process. We have identified t 0 roughly with the coda time of an
earthquake (Kagan & Knopoff, 1981), and have scaled this time as the square root of the
scalar seismic moment.
The exponent -3 /2 also arises in studies of Brownian motions. Thus, heuristically, we can
argue for a model of fracture as follows. At the moment that any given elementary earth-
quake rupture stops, the stress near the edge of the crack is lower than the critical breaking
stress for extension of the fracture. Let us assume the subsequent stress history near the
crack tip is due to a set of random factors due to other fractures in the neighbourhood. In
this case, the time history of the stress might resemble a Brownian random walk. The stress-
 Random stress and earthquake statistics 725
time function is thus given by the solution to the diffusion equation. When the stress reaches
the critical threshold level, a new earthquake begins. The distribution density of the time
intervals of 'first passage' (Feller 1966) is the function
a
f(t)= ·exp [-a 2 /(2Dt)], (I)
(2rrDt 3 ) 112
where D is the diffusion coefficient, t is the time, and a is the threshold or barrier stress, i.e.
the difference between the initial breaking stress and the stress at the time of cessation of
extension. This distribution as a function of stress is the Rayleigh distribution; as a function
of time it is the Levy distribution ( cf. Zolotarev 1983, p. 79). In this model the stress is
taken to be a scalar, which corresponds to the addition of perfectly aligned stress tensors;
the interaction of misaligned tensors will be discussed in a second contribution. The Levy
distribution (1) has been used by Kagan (1982) [see equations (8) and (9)] to model the
time dependent part of the space-time sequence of these infinitesimal events. It is only the
time sequence that concerns us here. The average recurrence time f is infinite for the model
of (1 ).
This model has been used to describe aftershocks, and in the aftershock case, as noted, it
yields the familiar Omori t- 1 aftershock frequency law. However, the usual time-scale of
aftershock series is of the order of months; this interval is too short to observe tectonic
effects. In fact, we can easily assume without penalty that the amount of stored energy
decreases monotonically during an aftershock sequence. On the other hand, the motion of the
plates restores the energy lost in the main shock and its subsequent aftershocks, in
anticipation of the next major earthquake. Thus, over a long time span, the statistics of
earthquakes must have a different interval-between-earthquakes law than is appropriate for
aftershocks. In this note, we discuss the appropriate pairwise interval law for a model in
which a steadily increasing source of stress, which we assume is due to plate tectonics, is
added to a random or diffusion component, where the distribution (1) describes the density
of earthquake recurrence times in the absence of tectonic loading. Under the model of
Brownian motions in the presence of an external steady field, the frequency of occurrence
law will be asymptotic to the Omori law for short time intervals. There is no guarantee that
the familiar exponential law for the time intervals between Poisson random events will be
the result.
If the rate of tectonic loading is a constant C, the distribution density f(t) is modified to
(Feller 1966,p.368)

a
f(t) = (2rr Dt3 )1/2 • exp
[--(a Wt
- Ct) 2
]
(2)

The average recurrence time in this case is no longer infinite, but becomes
1/2

f=~± ( ~~) (3)

where we have also given the standard deviation.

T/
If we set ~ = a/v'2JJ and =C/Vib we obtain the distribution density in terms of only
two parameters,

f(t) = Jrrt- ! 3 2
• exp - ( .Jr-T/Vt r (4)
726 Y. Y. Kagan and L. Knopoff

)-
....
~

~. 0100+-~~-+<f-H~~+-~~++-~~~+-+--+__.__.._....,~~f--~~t-+<-~~--1
w
0
)-
....
...J
~. 0010+-~~-Hf'~~~+-~-+---+~~~t--~~-tt-~-+-~t--~"'rl-+-t--~--1
<I

"'
0
""
IL

. 1 0 1. 0" 1 0. 0 0 100. "0

TIME

B
e--720°
u - c
71 -720 CRITICAL
--+ Tt.1E

~~ms~

..,,=o
Figure l. Time interval distributions. The probability densities for equations (1) and (4) are displayed for
various values of the dimensionless random ~ and external stress 11 fields. Time is measured in arbitrary
units which depend on the values of the quantities D and C. In the lower part of the figure we show
schematically Brownian motions in the absence and presence of an external stress field.

It is a trivial matter to rewrite (3) in these variables. If we pass to the limit 11---+ 0, we get ( 1)
which can be written as a function of~ and t.
The time dependence of the stress is governed by the Fokker-Planck equation (Feller
1966). Depending on the relative sizes of~ and 71, i.e. on the scaled value of the tectonic
loading rate, the resulting distribution may vary from the power-law or scale-invariant law
to a Gaussian renewal distribution (Fig. 1). The curve with ~ = 1, 11 = 0, for example,
corresponds to equation (1), and has a power-law long-time tail; as already noted, this is
appropriate for the case of aftershocks. As the value of the normalized drift velocity 11
increases, the resulting distribution of earthquake intervals becomes asymptotically
Gaussian. Similarly, if the normalized stress drop~ increases for non-zero 71, the distribution
of earthquake intervals also becomes almost Gaussian (see Fig. I).
The plots at the bottom of Fig. I illustrate schematically several realizations of the time
histories of randomly changing stresses for the cases of absence (left) and presence (right) of
tectonic loading. The stress trajectories start from a stress level which is less than the critical
threshold stress; under the influence of random factors it subsequently increases or
decreases. An earthquake starts at the time that it reaches the critical breaking stress.
 Random stress and earthquake statistics 727

0.001

0.00011---~~-+-~~-r-~~--+-~~--+~~~+---~

1 10 100 1000 t:l,000 -00,000

Time
Figure 2. Numbers of dependent earthquakes in several time intervals for the USGS catalogue of earth-
quakes for each of four (numbered) fixed distance intervals from the main shock (see equation 5 ).
Successive time intervals increase by a factor of .jiO.

We find partial confirmation of the above model from the time distribution of the
numbers of dependent shocks occurring in certain time and distance intervals from other
shocks (Kagan & Knopoff l 986a) for the USGS Central California catalogue (Fig. 2). The
distribution density µ(t) is the number of dependent earthquakes that occur in the
magnitude interval between M and M - 0.1, where M is the magnitude of a main shock. (In
practice, the number of aftershocks, or dependent shocks with magnitudes near Mis close to
zero. To determine this quantity we smooth the distribution of the number of aftershocks
as a function of magnitude to determine this asymptotic value. The smoothing makes use of
the function (I 0/(3) • µ · exp [(3(M - Ma)] analogous to the magnitude-frequency law, where
Ma is the magnitude of the dependent shock, and (3 is an adjustable parameter which is close
to In 10 for most catalogues.) Successive time intervals in Fig. 2 increase logarithmically; the
µ values are per unit logarithmic time interval, instead of per unit linear time interval as in
Fig. 1. The unit time in Fig. 2 is 2.22 x 10-4 days, which is our estimate of the coda length
for a main earthquake with ML= 1.5. Because of the logarithmic time scale, a flat curve in
Fig. 2 implies a power-law distribution. Since the roll-off in Fig. 2 is about 10 4 time units,
we infer that for times less than about 2 days, the temporal properties of aftershock
sequences are dominated by local random effects and for times greater than 2 days, tectonic
influences begin to be felt, for main earthquakes with ML = 1.5. The unit time is scaled
according to the magnitude of the main shock as
TM= 2.22x10-4 • 3M-l.5_

Details of smoothing and scaling procedures are to be found in Kagan & Knopoff ( 1980,
1986a).
The constancy of any of the curves in Fig. 2 indicates that the number of dependent
shocks falls off as l/t. However, this roughly uniform rate drops drastically for large time
intervals. We have used a fairly complicated spatial scaling relation because of the presence
728 Y. Y. Kagan and L. Knopoff
of location errors in the catalogue. The formula for the distance boundaries is
Ri = [r/ + (Pj • 3M-4.o )2 ]1/2, (5)

where ri = 0.35, 0.70, 1.1, 1.6 km, and Pj = 0.50, 1.25, 2.5, 5.0 km. (For the more spatially
distant curves there is a weak maximum that is shifted towards later times; this result
suggests an outward 'migration' of seismicity.) For comparison, the smooth curve of Fig. 2 is
drawn for the case of a normalized version of equation (4) with ~ = 0, 1/ =I 0- 2 - 25 • We
observe that a small but non-zero value of the drift component 1/ causes the power-law
distribution density of (1) to transform to an exponential distribution for large time
intervals and is similar to the results for real aftershock sequences as well. The time rate of
occurrence of aftershocks in other earthquake catalogues is similar (Kagan & Knopoff 1980).
For~}> T/, a good approximation to equation (4) is the gamma distribution (cf. Kagan &
Knopoff l 986a, section 2)
f(t)=O t.;;,t 0
(6)
2
where t x = 11- is a parameter that controls the distribution for large values of time and t 0 is
the characteristic dead time of the earlier model, which we identify as t 0 =4~ 2 /n. The
coefficient c is
2 2yn
c = 2 • t 01l 2 exp (-t 0 /tx) - 4 • r; 1! 2 Erfc [(t 0 /tx )1/ 2 ] t. -
'.::::'.. . t.,
vto vtx
where Erfc (x) is the complementary error-function. For time intervals close to t 0 , formula
( 6) is not a particularly good approximation to ( 4); for large time intervals the agreement
between the formulas is very good.
The above discussion is valid for sequences of the strongest earthquakes in a region. The
source volume of a weak earthquake will most likely be imbedded in the focal zone of a
future stronger event; thus the renewal chain for the specific volume will be destroyed. For
the strongest earthquakes, the stress levels in formulas (I)-( 4) are to be scaled according to
the size of the focal region. In effect, we interpret a in the above equations as an
appropriately scaled integral of the stress-drop over spatial extent of the earthquake focal
region. Thus a is rather close to the scalar moment of the earthquake; hence, the recurrence
time in (3) will be dependent on the size of event. The precise way to scale this relation is
not obvious. It is, of course, exactly this scaling that yields the seismic moment-frequency
law. We have discussed the time scaling of the seismic moment in Kagan & Knopoff ( 1981);
the spatial scaling of a focal region has been discussed by Kagan (1982). The distribution (2)
depends on two parameters, so its scaling should be more complicated than that for (I),
discussed in Kagan & Knop off ( 198 I).
As noted, another limitation of the discussion in this paper is the quasistatic nature of the
model; we do not consider here dynamical rupture effects which, most probably, control the
final size of an earthquake focal zone.
The model we have used is a rather crude approximation to the real picture. In reality,
stress is a tensor, rather than a scalar. The deviation of the seismic moment tensor from
complete coherence during the process of rupture may be the cause of intermittency of
rupture (Kagan & Knopoff 1985; we will discuss this point extensively in our companion
paper, see Kagan & Knopoff l 986b); most probably, these disorientations of the moment
tensor are due to asperities or barriers in the propagation of a fault. To study these problems
we need to consider the properties of the stress tensor, and the spatial distribution of the
elementary events.
 Random stress and earthquake statistics 729
3 Discussion

These demonstrations show that we can obtain distributions of interearthquake time

intervals that vary from a power-law to an almost Gaussian renewal type distribution. Both
of these limiting distributions are the result of rather simple assumptions concerning the
behaviour of the stresses in the lithosphere. The second of these limiting distributions is of
the 'repulsive type' and has been used by many advocates of the 'characteristic earthquake'
hypothesis (see, for instance, Wesnousky et al. 1983; Schwarz & Coppersmith 1984) to
model the distribution of main sequence earthquakes. Unfortunately, these repulsive
distributions of the interoccurrence time, which are not of cluster type, are not readily
adaptable to the multidimensional generalization which is necessary to describe spatial
distributions as well. It is no accident that all models based on repulsive-type renewal
distribution are 1-D. To reduce the dimensionality of the earthquake process to one time
dimension only, we need to isolate the region under consideration, a procedure which is
difficult to formalize. Unfortunately, opinions regarding methods of spatial isolation are
usually divergent. Put simply, we do not know how to minimize edge effects or the effects
of earthquakes in adjacent or other more distant regions, on the statistics of an extended
region considered as a point in space.
The (t ri) distributions cannot be directly compared to the known distributions of inter-
earthquake time intervals or to the distributions of numbers of aftershocks following a
main earthquake, because the clustering of earthquakes must be taken into account. ln
simulations of earthquake occurrence (Kagan & Knopoff 1981 ), we model this clustering
by a hypothesis that earthquakes are composed of infinitesimal events which multiply
according to a critical branching process. We show that the power-law distribution of inter-
event time with the exponent value of -3/2, combined with the above branching
assumption, reproduces aftershock-foreshock properties of shallow seismicity. These
properties are studied by an application of the same statistical analysis technique to both
real and synthetic catalogues (Kagan & Knopoff 1981).
We can also offer the following heuristic basis for the above agreement. As mentioned in
the introduction, the best known of the statistical distributions for temporal occurrence is
the Omori law of aftershock occurrence. This law is usually formulated as
N(t, t:.t) o: t-a • !:.t, (7)
where N(t) is the number of aftershocks in a small interval of time (t, t + !:.t), t is the elapsed
time after a main event, and a is an exponent close to one.
The distributions (1) and (2) have been calculated for interevent time intervals; we need
to compare them with the distribution of the numbers of earthquakes. We use the fractal
model of Mandelbrot (1982) to make the transition from the former distribution to the
latter.
If we assume that earthquakes form a temporal, 1-D, renewal chain of events, then our
results (equation 1) suggest that the exponent a should be 1/2; in other words the fractal
dimension of an earthquake sequence should be 1/2. The fractal dimension is defined for our
purposes as the exponent in the dependence of the cumulative number of events on the
time, which is one more than the exponent on the time in (1 ), for sufficiently large time-
intervals (Mandlebrot 1982). The critical branching process which we use to approximate
earthquake ruptures (Kagan & Knopoff 1981 ), has a stationary first moment; thus, in
principle, it can be regarded as a renewal chain of events. However, two considerations must
be taken into account before comparing our results with the Omori law: (1) higher moments
of the critical branching process are not stationary, and (2) the numbers in (7) are not
730 Y. Y. Kagan and L. Knopoff
counted from any arbitrary event, as suggested by Mandelbrot's (1982) definition of the
dimension, but instead are conditional on the occurrence of the especially strong cluster of
elementary events that comprise a main earthquake.
The total number of events at the beginning stages of development of a critical branching
process is proportional to n 2 (Athreya & Ney 1972, p. 20), where n is the number of
generations in the process. The number of new events is proportional ton. Thus, the fractal
dimension of the set of events at the beginning of the sequence should be ( l /2) x 2 = 1, (see
Mandelbrot 1982, eh. 32), i.e. there should be a continuous stream of events at the
beginning of any earthquake sequence. This stream of events is what we call an earthquake:
the exact definition of the time of termination of an earthquake depends on the frequency
response of the seismographic network, as well as other properties of the network (see more
in Kagan & Knopoff l 986a). New episodes of energetic branching in the same sequence
should be interpreted as new, separate earthquakes (aftershocks). When the sequence of
events starts to die out, the numbers of new events should first stabilize at an interval rate
-3/2, and then fall off to zero. Thus, the fractal dimension of the set should fall first to 1/2,
and then to zero. The zero value of the dimension corresponds to a= 1 in ( 7), i.e. the
cumulative number of events in the time interval T increases only as log (T).

4 Conclusions
The Omori law of aftershock occurrence and, in general, the time clustering of earthquake
events, has been shown to be easily derived from a simple assumption of the behaviour of
the stress prior to an earthquake and from elementary probability theory. If tectonic stress
loading is present, then the interearthquake time distribution may change from a power-
law to an almost Gaussian renewal distribution, depending on the ratio of the random
component to the drift component.

Acknowledgments
This research was supported in part by Grant CEE-84-07553 of the National Science
Foundation.
Publication Number 2893, Institute of Geophysics and Planetary Physics, University of
California, Los Angeles, CA 90024.

References
Athreya, K. B. & Ney, P. E., 1972. Branching Processes, Springer-Verlag, New York.
Feller, W., 1966. An Introduction to Probability Theory and its Applications, 2, J. Wiley, New York.
Kagan, Y. Y ., 1982. Stochastic model of earthquake fault geometry, Geophys. J. R. astr. Soc., 71,
659-691.
Kagan, Y. Y. & Knopoff, L., 1980. Dependence of seismicity on depth, Bull. seism. Soc. Am .. 70,
1811-1822.
Kagan, Y. Y. & Knopoff. L., 1981. Stochastic synthesis of earthquake catalogs, J. geophys. Res., 86,
2853-2862.
Kagan, Y. Y. & Knopoff, L., 1985. The two-point correlation function of the seismic moment tensor,
Geophys. J. R. astr. Soc., 83, 636-657.
Kagan, Y. Y. & Knopoff, L., l 986a. Statistical and stochastic models of earthquake occurrence and earth-
quake prediction (in preparation).
Kagan, Y. Y. & Knopoff, L., 1986b. Random stress and earthquake statistics: spatial dependence (in
preparation).
Mandelbrot, B. B., 1982. The Fractal Geometry of Nature, W. H. Freeman, San Francisco.
 Random stress and earthquake statistics 731
Schwarz, D. P. & Coppersmith, K. 1., 1984. Fault behaviour and characteristic earthquakes: Examples
from Wasatch and San Andreas fault zones, J. geophys. Res., 89, 5681-5698.
Wesnousky, S. G., Scholz, C. H., Shimazaki, K. & Matsuda, T., 1983. Earthquake frequency distribution
and the mechanics of faulting, J. geophys. Res., 87, 9331-9340.
Zolotarev, V. M., 1983. Odnomernye ustoichfrye raspredeleniya (One-Dimensional Stable Distributions),
Nauka, Moscow, (in Russian).