Geophys. J. Int.
(2010)
183,
1014–1030
doi: 10.1111/j.1365-246X.2010.04774.x
GJI Seismology
Probabilistic prediction of rupture length, slip and seismic ground
motions for an ongoing rupture: implications for early warning
for large earthquakes
Maren B
̈
ose and Thomas H. Heaton
Seismological Laboratory, California Institute of Technology (Caltech),
1200
E. California Blvd., Pasadena, CA
91125
, USA. E-mail: mboese@caltech.edu
Accepted 2010 August 12. Received 2010 August 11; in original form 2010 February 12
SUMMARY
Earthquake Early Warning (EEW) predicts future ground shaking based on presently available
data. Long ruptures present the best opportunities for EEW since many heavily shaken areas
are distant from the earthquake epicentre and may receive long warning times. Predicting
the shaking from large earthquakes, however, requires some estimate of the likelihood of the
future evolution of an ongoing rupture. An EEW system that anticipates future rupture using
the present magnitude (or rupture length) together with the Gutenberg-Richter frequency-
size statistics will likely never predict a large earthquake, because of the rare occurrence of
‘extreme events’. However, it seems reasonable to assume that large slip amplitudes increase
the probability for evolving into a large earthquake. To investigate the relationship between the
slip and the eventual size of an ongoing rupture, we simulate suites of 1-D rupture series from
stochastic models of spatially heterogeneous slip. We find that while large slip amplitudes
increase the probability for the continuation of a rupture and the possible evolution into a
‘Big One’, the recognition that rupture is occurring on a spatially smooth fault has an even
stronger effect. We conclude that an EEW system for large earthquakes needs some mechanism
for the rapid recognition of the causative fault (e.g., from real-time GPS measurements) and
consideration of its ‘smoothness’. An EEW system for large earthquakes on smooth faults,
such as the San Andreas Fault, could be implemented in two ways: the system could issue
a warning, whenever slip on the fault exceeds a few metres, because the probability for a
large earthquake is high and strong shaking is expected to occur in large areas around the
fault. A more sophisticated EEW system could use the present slip on the fault to estimate the
future slip evolution and final rupture dimensions, and (using this information) could provide
probabilistic predictions of seismic ground motions along the evolving rupture. The decision
on whether an EEW system should be realized in the first or in the second way (or in a
combination of both) is user-specific.
Key words:
Probabilistic forecasting; Earthquake dynamics; Earthquake ground motions;
Seismicity and tectonics; Statistical seismology; Early warning.
INTRODUCTION
Earthquake Early Warning (EEW) requires fast and robust predic-
tions of earthquake source and ground motion parameters shortly
after the initiation of an earthquake. If given in a timely manner,
warnings can be used to trigger and execute automatic measures to
reduce expected damage at distant sites (e.g. Goltz 2002; Allen
et al.
2009). Several algorithms for EEW have been proposed and tested
during recent years (e.g. Nakamura 1988; Allen & Kanamori 2003;
Kanamori 2005; B
̈
ose 2006; B
̈
ose
et al.
2007; Cua & Heaton 2007;
Wu
et al.
2007; Hoshiba
et al.
2008; B
̈
ose
et al.
2009a,b). The ma-
jority of these approaches use the amplitude and frequency content
of seismic waves at one or more seismic sensors in the epicentral
area for a rapid assessment of earthquake source and ground motion
parameters. With a few exceptions (e.g. Yamada
et al.
2007; B
̈
ose
et al.
2008; Yamada & Heaton 2008), the majority of the current
EEW approaches consider the earthquake as a point source, that is,
neglect rupture finiteness.
Large magnitude earthquakes (
M
>
7.0), with rupture lengths
of up to hundreds of kilometres, cause damaging ground shaking
in much larger areas than moderate-to-strong events. Despite their
rare occurrence, many more people and users could benefit from an
EEW system during large earthquakes (Heaton 1985; Allen 2006).
Due to the low propagation speed of seismic ruptures (
∼
80 per cent
of the seismic
S
-wave speed), warning times to heavily shaken areas
may exceed more than 1 min, while users who are strongly shaken
1014
C
2010 The Authors
Geophysical Journal International
C
2010 RAS
Geophysical Jour
nal Inter
national
Early warning for large earthquakes
1015
in smaller earthquakes are likely to be near the epicentre and will
generally receive warnings of a few seconds only.
A key to EEW for large magnitude earthquakes is the rapid
prediction of the potential final dimensions of the ongoing rupture.
This is because the level and distribution of seismic ground motions,
such as peak ground velocity (PGV) or seismic intensity, are usually
predicted from the earthquake magnitude and the rupture-to-site
distance, which both can be estimated from the expected rupture
dimensions. If the rupture is still propagating, these predictions
can only be given in a probabilistic manner. Appropriate decision
making is necessarily based on probabilistic estimates of future
shaking.
The rupture of a large (strike-slip) earthquake propagates mainly
in the horizontal direction along a fault, that is, rapid prediction of
the rupture length is especially important for EEW. An empirical
relationship between the subsurface rupture length
L
(km) and the
(moment) magnitude
M
for strike-slip earthquakes with
M
4.8 to
M
8.1, for example, has been proposed by Wells & Coppersmith
(1994) with
M
=
4
.
33
+
1
.
49 log(
L
)(1)
[see Working Group on California Earthquake Probabilities (1999)
for discussions on this relation].
A simple way to predict
L
(without any further information on the
rupture) is by the usage of frequency-size statistical distributions
of earthquakes, such as produced by the Gutenberg-and-Richter
model (G-R; Gutenberg & Richter 1944) or by the Characteristic
Earthquake model (CE; Schwartz & Coppersmith 1984).
The G-R model describes the frequency-size distribution of earth-
quakes by the power-law
log
N
(
M
)
=
a
−
bM
,
(2)
where
N
(
M
) is the number of earthquakes exceeding magnitude M,
and
a
and
b
are constants (Gutenberg & Richter 1944). The
b
-value
in eq. (2) is usually observed to be about 1.0, although there are
regional variations (e.g. Wiemer & Wyss 2002; Schorlemmer
et al.
2004);
b
also depends on the particular definition of magnitude
chosen (e.g.
M
l
,
M
S
and
M
w
).
The probability of occurrence of an earthquake exceeding a cer-
tain magnitude
M
in the range [
M
min
,
M
max
] can be described by
the probability density function (PDF)
prob
G
−
R
(
M
)
=
N
(
M
)
∫
M
max
M
p
N
(
M
)
d
M
=
β
e
−
β
M
e
−
β
M
min
−
e
−
β
M
max
=
β
e
−
β
(
M
−
M
min
)
1
−
e
−
β
(
M
max
−
M
min
)
,
(3)
with
β
=
b
ln[10] and
∫
M
max
M
min
prob
G
-
R
(
M
)d
M
=
1 (Fig. 1a). The
probability of exceedance (POE) of a present magnitude
M
p
is
p
G
−
R
(
M
|
M
p
)
=
prob
G
-
R
(
M
)
prob
G
-
R
(
M
p
)
=
e
−
β
(
M
−
M
p
)
=
10
−
b
(
M
−
M
p
)
.
(4)
If, for example, we knew that the present magnitude of an ongoing
earthquake was
M
p
5.5 (corresponding to a present rupture length of
L
p
≈
6 km) and
b
=
1, then the G-R model predicts that there is only
a 1 per cent chance that the event will grow to be larger than
M
7.5
(Fig. 1b). An EEW system based on the G-R model will always
conclude that there is a very low probability that such an earthquake
will grow to be a ‘Big One’; only 1 out of 100 earthquakes with a
present magnitude of
M
p
5.5 will eventually grow to be larger than
M
7.5. Assuming higher
b
-values in the G-R model will increase
the probabilities of earthquake growth (Iwata
et al.
2005).
The CE model, in contrast, assumes that single faults tend to
produce earthquakes of similar size. The magnitudes of these ‘char-
acteristic’ earthquakes mainly depend on the dimension of the fault
zones or fault segments (Schwartz & Coppersmith 1994). It is typ-
ically supposed that characteristic earthquakes are accompanied by
foreshocks, aftershocks and low-level background seismicity along
the considered fault or fault zone, which are distributed according to
the G-R model. Youngs & Coppersmith (1985) and Convertito
et al.
(2006) give a PDF of the earthquake occurrence that is appropriate
for this model with
prob
CE
(
M
)
=
⎧
⎨
⎩
β
e
−
β
(
M
−
M
min
)
1
−
e
−
β
(
M
max
−
M
min
−
M
2
)
1
1
+
c
for
M
min
≤
M
<
M
c
β
e
−
β
(
M
max
−
M
min
−
M
1
−
M
2
)
1
−
e
−
β
(
M
−
M
min
−
M
2
)
1
1
+
c
for
M
c
≤
M
≤
M
max
,
(5)
where the constant
c
is given by
c
=
M
2
β
e
−
(
M
max
−
M
min
−
M
1
−
M
2
)
1
−
e
−
β
(
M
max
−
M
min
−
M
2
)
.
(6)
Figure 1.
(a) Probability density function (PDF) of magnitude
M
and (b) probability of exceedance (POE) of different
M
for a rupture with a present magnitude
of
M
p
5.5. The probability for the occurrence of a large earthquake is much lower in the Gutenberg-and-Richter (G-R) than in the Characteristic Earthquake
(CE) model: only 1 out of 100 earthquakes will exceed
M
7.5. An EEW system based on the G-R model will likely never predict a ‘Big One’. The CE model,
on the other hand, is very poorly constrained (here:
b
=
1.0,
M
1
=
1.75 and
M
2
=
0.5).
C
2010 The Authors,
GJI
,
183,
1014–1030
Geophysical Journal International
C
2010 RAS
1016
M. B
̈
ose and T. H. Heaton
The parameters
M
1
and
M
2
describe two intervals below and
above the characteristic magnitude
M
c
=
M
max
–
M
2
(Fig. 1a).
The POE for the CE model,
p
CE
, can be determined analogous
to eq. (4). Due to the occurrence of characteristic earthquakes,
p
CE
takes generally higher values than
p
G
-
R
. For example with
M
1
=
1.75 and
M
2
=
0.5, the CE model predicts that there is a 40 per
cent chance for an ongoing event with
M
p
5.5 to grow to be larger
than
M
7.5, compared to 1 per cent predicted by the G-R model
(Fig. 1b).
Spatially smooth ruptures along so-called ‘mature’ faults usu-
ally travel over larger distances than heterogeneous ruptures along
‘immature’ or ‘generic’ faults (e.g. Kanamori & Allen 1986;
Wesnousky 1988; Manighetti
et al.
2005, 2009). For this study
we assume that the frequency-size statistics of earthquakes along
mature faults can be described by the CE model, while the G-R
model seems to be more appropriate to characterize the seismicity
along generic faults. There is an ongoing debate among seismolo-
gists about which model might be more appropriate to characterize
the seismic activity along individual faults (e.g. Wesnousky 1994;
Kagan 1996; Page
et al.
in preparation). We want to avoid these
discussions here. Even if we knew that the G-R or CE models are
adequate to describe the seismicity in a certain region or along a
specific fault, it seems clear that there is insufficient knowledge to
assign appropriate parameters to the wide range of possible rupture
scenarios.
A more flexible approach to study the rupture probabilities along
faults is by the application of stochastic slip models that can be
tuned from G-R- to CE-like behaviour. These models can be used
to generate series of spatially variable slip (Liu-Zeng
et al.
2005).
Slip heterogeneities can be caused, for example, by the frictional
properties and the spatial heterogeneities in the material character-
istics on the fault, or by the complexities of the fault geometry (Rice
1993; Andrews 1994; Aki 1995; Aagaard & Heaton 2008). In this
study, we use 1-D stochastic slip models to simulate large suites of
possible rupture scenarios for a systematic statistical analysis of the
relationship between seismic slip and rupture length.
We find that the larger the present slip amplitude
D
p
, the higher
is the probability for the continuation of a rupture and the evolution
into a ‘Big One’, that is, a long remaining rupture length
L
r
(Yamada
2007). We also find that the
D
p
–
L
r
relationship shows a significant
amount of scatter and is strongly controlled by the degree of assumed
spatial heterogeneity of slip on the fault.
In this paper, we are not addressing the origin of spatially het-
erogeneous slip. We simply note that if slip changes rapidly, then
it is more difficult to predict the future evolution of the rupture (in
space and time) than if the slip is spatially smooth. In principle,
we could also try to estimate the smoothness of an ongoing rupture
in real time from the variability of slip amplitudes relative to the
mean evolution of slip. This idea will be briefly discussed in the
next section, but shall not be a major topic in this study.
The real-time application of our probabilistic approach to EEW
requires the continuous monitoring of seismic slip along faults, for
example, by a global positioning system (GPS; e.g. Crowell
et al.
2009) or by the backprojection of seismic displacement data onto the
fault line (Yamada 2007). Again, to understand the implications of
a small or large slip amplitude, it is essential to consider the respec-
tive ‘smoothness’ of the fault along which a detected earthquake
propagates. This means that an EEW system for large earthquakes
needs some mechanism for the rapid recognition of the rupturing
fault and the capability to consider the respective statistical features.
This procedure requires that the rupture is automatically linked to
a (hopefully) known geological structure, for which a ‘smoothness
measure’ from an earlier assessment was determined and stored in
a database.
If rupture is occurring on a smooth fault, for example, on the San
Andreas Fault (SAF), an EEW system could be implemented in two
ways: the system could issue a warning, whenever slip on the fault
exceeds a few metres, because the probability for a large earthquake
is high and strong shaking is expected in large areas around the
fault. A more sophisticated EEW system could use the present slip
on the fault to estimate the final rupture dimensions and future slip
evolution, and (using this information) could provide probabilistic
predictions of seismic ground motions along the evolving rupture.
We will demonstrate this concept of an EEW system on the example
of the southern Californian SAF.
DATA
1-D stochastic slip models
We assume that seismic ruptures consist of slip pulses that propagate
along a fault (Heaton 1990). The ultimate slip
D
(
x
) at the spatial
point
x
on the fault is achieved quickly after the passage of the
rupture front. That is, we suppose that the current distribution of
slip up to a present rupture length
L
p
is approximately the final
(static) distribution of slip along this length, and no further slip
accumulation occurs after the passage of the rupture front.
Although earthquakes actually occur on rupture surfaces, we are
constructing here a simple simulation of long ruptures and assume
a 1-D line source. We assume that slip pulses vary (perhaps chaot-
ically) as they propagate along the fault. We also assume that the
event is over, once the amplitude of the slip pulse drops to zero.
We simulate the propagation of slip pulses from simple 1-D
stochastic slip models as proposed by Liu-Zeng
et al.
(2005), in
which the slip amplitudes at adjacent points,
x
and (
x
+
x
), change
by some random amount in either positive or negative direction. We
model this change in the amplitude by a low-pass filtered Gaussian
white noise series
R
(
x
) with zero mean and standard deviation one,
D
(
x
)
=
D
|
FT
−
1
[
ˆ
R
(
k
)
k
−
μ
]
|
,
(7)
where
ˆ
R
(
k
) is the Fourier transform of
R
(
x
),
k
is the wavenumber,
FT
−
1
the inverse Fourier transform and
D
a scaling factor. The
filtering exponent
μ
determines the spatial roughness of slip: the
smaller
μ
, the rougher is the slip. Each point
x
,atwhich
D
(
x
)
=
0, defines the end of an individual rupture. We take the absolute
amplitudes of these ‘parent series’ to obtain contiguous positive slip
functions (Fig. 2). Note that the slip models described by eq. (7) are
purely geometric; they only assume that slip is spatially contiguous
and that it has the scale-invariant properties of a self-affine fractal
(Turcotte 1997).
Liu-Zeng
et al.
(2005) found that models with a rupture roughness
of up to
μ
=
1.5 produce slip series that have similar statistical
features like real earthquakes (see below for further discussions).
We use the stochastic slip models in eq. (7) to generate a set of 500
parent series of 65 536 points each for
μ
=
1.1,
μ
=
1.3 and
μ
=
1.5. A total of 22 541 individual slip events, which are identified
by determining adjacent zero crossings of the parent series, are
catalogued: 15 006 for
μ
=
1.1, 5121 for
μ
=
1.3 and 2414 for
μ
=
1.5.
Models with small rupture roughness
μ
, that is, high slip hetero-
geneity, tend to generate more small events, that is, shorter ruptures,
than models with large
μ
(Fig. 3a), because rougher slip distribu-
tions are more likely to terminate rupture (Fig. 2). The frequency-
size statistics of these events are close to the power-law distributions
C
2010 The Authors,
GJI
,
183,
1014–1030
Geophysical Journal International
C
2010 RAS
Early warning for large earthquakes
1017
Figure 2.
Unscaled parent series of slip generated from eq. (7) for three
levels of rupture roughness
μ
. The larger the slip heterogeneity, the larger is
the probability to cross the zero line, which defines the end of an individual
rupture event.
generated by the G-R model; slip series generated from models with
large
μ
, on the other hand, produce smooth ruptures with a higher
number of large events, that is, the frequency-size statistics of these
series are more similar to the distributions produced by the CE
model (Fig. 3a). This means that we can use the roughness parame-
ter
μ
to tune the stochastic models from G-R- to CE-like behaviour.
This is in good agreement with, for example, Hillers
et al.
(2007),
who found from quasi-dynamic rupture simulations that smooth
continuum faults usually tend to have narrow event distributions,
and that strong spatial heterogeneity is required to produce spatio-
temporal non-uniform seismic behaviour.
The calibration of the 1-D slip functions is driven by two desires:
(1) we want to reproduce realistic average slip-to-rupture length
ratios as observed during earthquakes, and (2) we want to avoid
maximum slip amplitudes that seem unrealistic high (although we
are not really sure what is unrealistic). In this study, we set in eq. (7)
D
=
3 m. We assume that the 1-D slip amplitudes correspond to
the average slip taken over the fault width.
The slip amplitudes
D
in our models follow one-sided Gaussian
distributions (Fig. 3b) as a consequence of the assumed Gaussian
distribution of noise amplitudes in eq. (7). If desired, we could re-
place the Gaussian by any other distribution, for example, by the
logarithmic distribution. The slip amplitudes in the stochastic mod-
els increase with increasing rupture roughness, that is, decreasing
μ
(Fig. 3b). The largest slip amplitudes in our models are less than
30 m for
μ
=
1.3 and
μ
=
1.5, and, in a very few cases, grow
up to 65 m for
μ
=
1.1. The range of simulated slip amplitudes
seems to be fairly realistic if compared, for example, with the slip
models in the Finite-Source Rupture Model Database (SRCMOD
version 7, http://www.seismo.ethz.ch/srcmod/). Note, however, that
the inverted slip distributions in the SRCMOD are usually strongly
smoothed and the SRCMOD does not contain a sufficient number
of large events to allow for a statistical analysis as aimed for in our
studies.
We generally find a good agreement between the ratios of the
average slip
D
(taken over the entire rupture length
L
) to the rup-
ture length
L
in the simulated slip functions and earthquakes which
usually have
D
/
L
values within the wedge-shaped area in Fig. 4:
the
D
/
L
ratios, which are related to the seismic stress drops, in-
crease among small earthquakes compared to large events, if we
consider the combined effect of different levels of rupture rough-
ness (Fig. 4). Ruptures with rougher slip distributions (smaller
μ
)
tend to produce higher slip-length ratios than smooth ruptures. Gen-
erally, the ratios decrease with increasing rupture length
L
,ex-
cept for smooth ruptures (
μ
=
1.5) that hardly show any variation
in
D
/
L
.
For rupture lengths that are large relative to the seismogenic
width of the fault (
L
30 km), constant slip versus length scal-
ing (of smooth ruptures) implies constant stress drops with
L
.This
behaviour is compatible with the ‘
L
-model’ proposed by Scholz
(1982). Note, however, that the stress drop is not really a fundamen-
tal frictional parameter in our 1-D stochastic slip models: the end of
a rupture is simply defined by the location at which the amplitude
of a slip pulse dies to zero. That is, the stress drop in our models is
determined by the details of spatially variable slip pulses (Liu-Zeng
et al.
2005).
Figure 3.
(a) Frequency-length and (b) frequency-slip statistics for the simulated slip events for three levels of rupture roughness
μ
using the 1-D stochastic slip
models proposed by Liu-Zeng
et al.
(2005). Smooth ruptures (larger
μ
) tend to generate more long ruptures and smaller slip amplitudes than rough ruptures.
The roughness parameter
μ
allows tuning the slip models from G-R- to CE-like behaviour. The range of simulated slip amplitudes agrees fairly well with the
slip amplitudes in the Finite-Source Rupture Model Database (SRCMOD version 7; http://www.seismo.ethz.ch/srcmod/) (dashed line). Note, however
,thatthe
inverted slip distributions in the SRCMOD are usually strongly smoothed and the SRCMOD does not contain a sufficient number of large events to allow for
a statistical analysis as anticipated in our study.
C
2010 The Authors,
GJI
,
183,
1014–1030
Geophysical Journal International
C
2010 RAS
1018
M. B
̈
ose and T. H. Heaton
Figure 4.
Average slip-to-rupture length
D
/
L
-ratios versus
L
obtained from
the 1-D stochastic slip models for different levels of rupture roughness
μ
.
The corresponding values of real earthquakes are typically located in the
wedge-edged area between the two lines. The scaling of the parent series
in eq. (7) is driven by the compromise between avoiding too large slip
amplitudes and the desire to simulate ruptures with fairly realistic
D
/
L
ratios.
Generally, the
D
/
L
ratios obtained from the simulated slip series
tend to be smaller than the ratios observed during real earthquakes
(Fig. 4, wedge-edged area). In principle, we could improve the fit
by increasing the scaling factor
D
in eq. (7), but this would lead
to single slip amplitudes exceeding several tens of metres, which
seems too large for a real earthquake. The decision to set in this
study
D
=
3 m is driven by the compromise between avoiding
slip amplitudes that are too large and the desire to simulate ruptures
with realistic
D
/
L
ratios.
Another way to evaluate the suitability of the simulated slip func-
tions is by comparing their shape with observations made during
real earthquakes. The mean shape of the 1-D slip functions can be
characterized by three parameters (Ward 2004): the (total) rupture
length
L
, the peak average slip
E
peak
and a skewness parameter
q,
E
(
D
L
,
E
peak
,
q
(
l
))
=
2
E
peak
l
q
L
2
q
(
L
2
q
−
l
2
q
)
1
/
2
,
(8)
where
q
controls whether a slip function has more off-centre (
q
≈
2.0) or more elliptical shape (
q
≈
0.5; Ward 1997). We determine
for each simulated slip series of length
L
parameters
E
peak
and
q
from fitting the model in eq. (8) to the generated 1-D slip functions
(Figs 5a–d).
The majority of the slip series tend to have elliptical shapes,
independent from the level of slip heterogeneity (Fig. 5c). This
symmetry is not surprising since the probability for slip increase
or decrease in the stochastic slip models is constant at all rupture
lengths
l
, that is, on average we should obtain the same shape,
independent from whether the rupture starts at
l
=
0orat
l
=
L
. Manighetti
et al.
(2005) found that the slip functions of real
earthquakes usually tend to be roughly triangular both along strike
and dip, and most of them (70–80 per cent) are asymmetric. The
reason for this phenomenon is not well understood, and cannot be
reproduced by our simple static 1-D models.
We have seen that the 1-D stochastic slip models as used in this
study imitate important statistical features of earthquakes, such as
reproducing realistic slip amplitudes and average slip versus rupture
length ratios (Fig. 4). On the other hand, the models fail to reproduce
asymmetric shapes of the slip series as often are observed during
real earthquakes. We think that employing the 1-D models for our
purposes, however, is justified, because we are mainly interested in
statistical features such as fractal behaviour of slip, which is well
reproduced by the 1-D functions (see Discussions and Conclusions).
What is the relationship between the simulated slip series
D
in
eq. (7) and the mean slip function
E
in eq. (8)? Neither the shape
factor
q
nor the peak slip value
E
peak
show a clear dependency on
μ
(Figs 5c and d). However, the variability of the simulated slip ampli-
tudes
D
relative to
E
increases strongly with increasing roughness of
the slip series, that is, the stronger the slip heterogeneity, the larger
are the differences of single slip amplitudes and the mean slip func-
tion, and the larger the standard deviation
σ
of the residuals (
E
-
D
)
(Fig. 6).
From our stochastic 1-D models we find that
σ
≤
1.5 m for
μ
=
1.5, and
σ
≤
2.25 m for
μ
=
1.1 and
μ
=
1.3 (Fig. 6). The variability
increases with increasing mean slip
E
. We could in principle use
the variability of the slip amplitudes along an ongoing rupture to
estimate the rupture roughness in real-time. This procedure, how-
ever, requires densely sampled slip profiles. In the following, we
assume that
μ
is a fault-characteristic parameter, which has been
determined from previous geological and geophysical studies (e.g.,
Sagy
et al.
2007) and is known when the rupture nucleates.
METHOD
Bayesian approach
The statistical analyses of seismic slip amplitudes
D
and rupture
length
L
in this paper are based on the Bayesian theorem. Proba-
bilistic methods, such as provided by the Bayesian framework, are
reasonable approaches to deal with the large variability in the ob-
served seismic source and ground motion characteristics and the
inherent uncertainties in predicting these parameters.
In this study, we determine the probabilities of remaining rup-
ture length
L
r
of an ongoing earthquake rupture using a presently
observed seismic slip amplitude
D
p
. We are not only interested in
determining the most probable, but also the less likely, however,
still possible solutions. In other words: we want to determine the
entire PDF for
L
r
. Both slip amplitudes
D
and rupture length
L
span several orders of magnitudes. They usually follow power-law
distributions, and we therefore use their logarithmic values log(
D
)
and log(
L
) in the following analyses.
We describe the rupture length–slip relationship by the ‘condi-
tional probability’ of
L
r
for a given
D
p
, also known as the ‘posterior
probability’
p
(log(
L
r
)
|
log(
D
p
)). The Bayesian theorem states that
the posterior probability can be determined from the product of the
likelihood function
p
(log(
D
p
)
|
log(
L
r
)) and the
apriori
probability
p
(log(
L
r
)), that is,
p
(log(
L
r
)
|
log(
D
p
))
=
p
(log(
D
p
)
|
log(
L
r
))
p
(log(
L
r
))
p
(log(
D
p
))
=
likelihood
×
apriori
normalization factor
,
(9)
where the normalization factor in the denominator,
p
(log(
D
p
))
=
n
i
=
1
p
(log(
D
p
)
i
|
log(
L
r
))
p
(log(
L
r
)), ensures
that the sum over all probabilities gives 1.
The likelihood function
p
(log(
D
p
)
|
log(
L
r
)) in eq. (9) describes,
how likely it is to observe a certain log(
D
p
)–log(
L
r
) combination.
C
2010 The Authors,
GJI
,
183,
1014–1030
Geophysical Journal International
C
2010 RAS
Early warning for large earthquakes
1019
Figure 5.
Statistical analyses of the shape of the simulated 1-D slip functions for three levels of rupture roughness
μ
. (a) and (b): Example slip series and mean
slip functions
E
(dashed lines);
E
depends on the peak mean value
E
peak
, shape factor
q
and rupture length
L
(see eq. (8)). (c) and (d): Relative frequency of
parameters
q
and
E
peak
for all 22 541 simulated ruptures. The majority of slip functions has a symmetric shape (
q
=
0.5) as a consequence of equal probabilities
of slip increase and decrease at all length
l
, that is, in a statistical sense the slip series should not be different if the rupture starts at
l
=
0or
l
=
L
.
Figure 6.
Standard deviation
σ
of residuals (mean slip
E
–slip
D
)forthree
levels of rupture roughness
μ
. The smaller the slip heterogeneity, the smaller
is
σ
, that is, the better is the fit between
E
and
D
.
The
apriori
probability
p
(log(
L
r
)) tells, how often we will observe
agivenlog(
L
r
) value. Obviously,
p
(log(
L
r
)) depends on both the
probability of occurrence of a rupture with length
L
andonhowfar
the considered rupture already has propagated, that is,
p
(log(
L
r
))
=
p
(log(
L
–
L
p
)), where
L
p
is the present rupture length, at which slip
D
p
is observed. In principle, we can use the log–linear relationship
between the rupture length and magnitude in eq. (1) and describe
the
apriori
probability
p
(
M
) by the earlier discussed frequency-size
statistics of earthquakes as produced, for example, by the G-R or
the CE model.
The likelihood function
p
(log(
D
p
)
|
log(
L
r
)) in eq. (9) can be deter-
mined from the joint probability function
p
(log(
D
p
),log(
L
r
)), which
describes the probability of the joint occurrence of log(
D
p
)and
log(
L
r
). The likelihood function
p
(log(
D
p
)
|
log(
L
r
)) can be calcu-
lated from
p
(log(
D
p
)
|
log(
L
r
))
=
p
(log(
D
p
)
,
log(
L
r
))
p
(log(
L
r
))
,
(10)
where we determine the joint probabilities from the 2-D kernel-
based estimation method (e.g. Bishop 1995) using Gaussian kernels
C
2010 The Authors,
GJI
,
183,
1014–1030
Geophysical Journal International
C
2010 RAS
1020
M. B
̈
ose and T. H. Heaton
with
p
(log(
D
p
)
,
log(
L
r
))
=
1
N
1
σ
√
2
π
×
N
n
=
1
exp
−
(log(
D
p
)
,
log(
L
r
))
−
(log(
D
p
)
,
log(
L
r
))
n
2
2
σ
2
,
(11)
where the standard deviation
σ
in the denominator controls the
width of the kernels, and
N
is the number of data points. Kernel-
based methods belong to the group of ‘non-parametric’ density
estimation techniques that do not require assumptions about the
forms of the PDFs (except that they are smooth). They therewith
are highly flexible and well suited for many applications (e.g. Bishop
1995).
RESULTS
Probabilistic prediction of rupture length
L
and magnitude
M
We go along each stochastically simulated 1-D slip series and de-
termine for each slip amplitude
D
p
the PDF of remaining rupture
length
L
r
, that is, the posterior probability
p
(log(
L
r
)
|
log(
D
p
)). We
consider two levels of slip heterogeneity in the stochastic models:
μ
=
1.1 (rough rupture) and
μ
=
1.5 (smooth rupture) represent-
ing ‘generic’ and ‘mature’ fault conditions, respectively. Note that
the slip and length values as specified in the following should be
understood in a qualitative rather than a quantitative way, because
they depend on the calibration of the 1-D slip models as described
in the previous section.
Fig. 7 shows the probabilities of remaining rupture length
L
r
for
slip amplitudes of
D
p
=
2.0 m,
D
p
=
5.0 m and continuous values
of
D
p
≤
20.0 m. The dashed lines mark the median values of the
PDFs, that is, 50 per cent of the simulated ruptures propagate over a
shorter and 50 per cent over a larger distance than specified by this
value. Note the strong asymmetry of the PDFs, leading to significant
differences between the median values and the maxima of the PDFs.
The greyish colours mark the 38 per cent, 68 per cent and 95 per
cent confidence intervals of the PDFs, corresponding to 0.5
σ
,1
σ
and 2
σ
of a Gaussian distribution.
Our model estimates a 50 per cent chance that a rough rupture
(
μ
=
1.1) keeps growing for
L
r
>
6kmif
D
p
=
2.0 m, and for
L
r
>
20 km if
D
p
=
5.0 m (Fig. 7, left-hand side). Using eq. (1)
these lengths correspond to magnitudes of
M
>
5and
M
>
6,
respectively. Note that a more accurate specification of
M
requires
the knowledge of the present rupture length
L
p
. The probabilities
for rupture continuation increase substantially, if the rupture travels
along a smooth (mature) fault with little slip heterogeneity (
μ
=
1.5;
Fig. 7, right-hand side): 50 per cent of the ruptures keep growing
for
L
r
>
45 km (
M
∼
7.0) if
D
p
=
2.0 m, and for
L
r
>
70 km (
M
>
7) if
D
p
=
5.0 m. Our model thus predicts that the potential sizes
of ruptures on heterogeneous (generic) and smooth (mature) faults
may differ by more than one full magnitude unit for the same slip
amplitude.
On a rough fault,
L
r
scales almost linearly with
D
p
(Fig. 7c, left-
hand side), while the confidence intervals become broader, that is,
the uncertainties in the prediction increase with increasing slip. Note
that even for large slip amplitudes (e.g.
D
p
=
6.0 m) every second
rupture does not exceed 25 km and therewith
M
<
6.5 (provided
that
D
p
is observed at
L
p
≈
0 km, i.e. close to the initiation point of
the rupture); for 38 per cent of the earthquakes ruptures we observe
L
r
<
50 km, for 68 per cent
L
r
<
75 km and for 95 per cent
L
r
<
125 km. If the rupture, however, occurs along a smooth (mature)
fault (Fig. 7c, right-hand side), we observe a strong increase of
L
r
for
D
p
≤
5.0 m. The PDFs hardly change for larger
D
p
amplitudes,
that is, the remaining rupture length is largely independent from the
present slip and every second rupture will grow for
L
r
>
65 km,
that is,
M
>
7.0.
Fig. 8 shows a comparison of the predicted probabilities of
remaining rupture length
L
r
and the relative frequency of
L
r
determined from the slip models in the Finite-Source Rupture
Model Database (SRCMOD version 7, http://www.seismo.ethz.ch/
srcmod/) for slip amplitudes of
D
p
=
1–5 m. Note that these plots
show observational slip data from both heterogeneous and smooth
ruptures. There is qualitative agreement between these data and the
probabilities derived from the stochastic models with slip hetero-
geneities of
μ
=
1.1 and
μ
=
1.3. The stochastic models capture
essential features of the SRCMOD. A more detailed comparison,
however, required a classification of the SRCMOD ruptures with
respect to their levels of slip heterogeneity, which needs to be done
in the future.
Prediction of mean slip function
Although the prediction of single slip amplitudes
D
along the evolv-
ing rupture may likely remain an elusive goal, we can try to estimate
the mean slip function
E
using eq. (8). This involves some knowl-
edge or the prediction of the total rupture length
L
, the peak mean
slip
E
peak
and shape factor
q
. The estimation of these parameters
requires some information on the evolution of slip. For example, if
we knew the entire slip series up to rupture length
L
p
, we could esti-
mate
L
=
L
p
+
L
r
from eq. (9) and then determine
E
peak
and
q
from
curve fitting of function
E
in eq. (8) to the observed slip amplitudes
D
. The shape of the slip function, however, remains elusive, if
E
peak
and
q
cannot be reasonably well resolved.
We test the robustness of the curve-fitting procedure to deter-
mine
E
peak
and
q
if the number of slip amplitudes is reduced by
downsampling of data points in the simulated slip series. We use
the slip amplitudes along the full rupture, that is, from
l
=
0to
l
=
L
, and assume that the total rupture length
L
of each event is known.
Fig. 9 shows the rms errors for both parameters as a function of data
completeness
C
. A completeness of
C
=
0.1, for example, means
that only every 10th slip amplitude of the original slip series is
used.
We can easily recognize a rapid decrease of the rms errors in
Fig. 9 with increasing data completeness
C
, that is, the more details
(slip amplitudes) of the slip series are known, the more reliable
are the estimates of
q
and
E
peak
. Rough ruptures (with small
μ
)
that have a large variability in the slip amplitudes require more
information on the evolution of slip than smooth series. Ruptures
with slowly varying slip maintain their memory of past rupture over
relatively longer lengths, that is, slip variations occur over larger
length scales. We call this behaviour the (spatial) ‘memory effect’
of smooth ruptures.
In the previous example we assumed that the rupture is already
terminated and slip amplitudes (or subsets of amplitudes) along this
rupture as well as the total rupture length
L
are known, before
q
and
E
peak
are determined from curve fitting. In the case of an ongoing
rupture, we will know the slip amplitudes only of the initial portion
of the rupture (if at all). The results in Fig. 9, however, indicate
that due to the (spatial) ‘memory effect’ of smooth ruptures the
prediction of future slip does not require detailed information on
C
2010 The Authors,
GJI
,
183,
1014–1030
Geophysical Journal International
C
2010 RAS
Early warning for large earthquakes
1021
Figure 7.
Probabilities of remaining rupture length
L
r
for a rupture with a present slip amplitude of (a)
D
p
=
2.0 m, (b)
D
p
=
5.0 m and (c) continuous values
of
D
p
≤
20.0 m. The probabilities are derived from the similated 1-D rupture series. Shown are the median values of the probability density functions (PDFs),
as
well as the 38 per cent, 68 per cent and 95 per cent confidence intervals, corresponding to 0.5
σ
,1
σ
and 2
σ
of a Gaussian distribution. Note that the probabilities
are plotted on a linear length scale, while determined on a logarithmic scale using eq. (9). Rough ruptures with a high level of slip heterogeneity (lef
t-hand
side) tend to propagate over shorter distances than smooth ruptures (right-hand side). The median values of the PDFs increase almost linearly with in
creasing
slip
D
p
for a rough rupture (c, left-hand side), while staying almost constant after a rapid increase up to
D
p
=
5.0 m for a smooth rupture (c, right-hand side).
the evolution of slip if a smooth (mature) fault is considered. We will
come back to this when demonstrating the probabilistic prediction
for a scenario earthquake along the SAF.
Implications for EEW
EEW requires a rapid and robust determination of seismic source
and ground motion parameters, such as magnitude and PGV, shortly
after the initiation of a moderate-to-large earthquake. Most chal-
lenging in EEW for large earthquakes is the real-time determination
of the probability that the rupture of an ongoing earthquake with a
present length
L
p
(corresponding to a present magnitude
M
p
) will
continue and exceed a certain magnitude level. We have seen earlier,
that in the G-R model the POE for a large magnitude earthquake is
extremely small: only 1 out of 100 earthquakes with
M
p
5.5 (
L
p
≈
6 km) will exceed
M
7.5 (Fig. 1). The output of the CE model, on
the other hand, depends very much on the model parameters chosen
and thus is not well constrained.
C
2010 The Authors,
GJI
,
183,
1014–1030
Geophysical Journal International
C
2010 RAS
1022
M. B
̈
ose and T. H. Heaton
Figure 8.
Comparison of the predicted probabilities of remaining rupture length
L
r
for two levels of slip heterogeneity
μ
with the relative frequency of
L
r
determined from the slip models in the Finite-Source Rupture Model Database (SRCMOD, http://www.seismo.ethz.ch/srcmod/); slip amplitudes vary f
rom
D
p
=
1–5 m from left to right. Note that these plots show observational slip data from both heterogeneous and smooth ruptures. There is good agreement
between these data and the probabilities derived from the stochastic slip models with
μ
=
1.1 and 1.3. A more detailed comparison required a classification of
the SRCMOD ruptures with respect to their levels of slip heterogeneity.
Figure 9.
Robustness of the curve-fitting procedure to determine (a) the shape parameter
q
and (b) the peak average slip
E
peak
of the mean slip function
E
in
eq. (8) if only a subset of slip amplitudes along the rupture is used. The rms errors decrease rapidly with increasing data completeness up to
C
=
0.2, which
means that every fifth slip amplitude is assumed to be known. The larger the slip heterogeneity, the larger are the errors, that is, the more important is
the
knowledge of past slip for predicting future amplitudes. Ruptures with slowly varying slip (higher
μ
) maintain their memory of past rupture over relatively
longer lengths. This is the so-called (spatial) ‘memory effect’ of smooth ruptures.
For the stochastic slip models the POE of a given magnitude
tends to increase with increasing slip amplitude
D
p
(Fig. 10). For
example, for
M
p
5.5 and a slip amplitude of
D
p
=
0.5 m it is three
times more likely that the rupture on a rough fault (
μ
=
1.1) will
become
M
>
7.0 and 10 times more likely to become
M
>
7.5
than in the G-R model (Table 1). The probabilities increase by an
additional factor of two, if the slip increases to
D
p
=
5.0 m. A
further increase of the probabilities by an additional factor of two
to three is observed, if the same slip occurs along a smooth fault
(
μ
=
1.5). This observation demonstrates that although the present
slip amplitudes
D
p
affect the POE, the smoothness of the causative
fault has an even larger impact.
When should an EEW system issue a warning that an earthquake
may become large? We have seen in Fig. 7c that
L
r
scales almost
linearly with
D
p
on a rough rupture (left-hand side) and even for
large slip amplitudes (e.g.
D
p
=
6.0 m) only every second rupture
exceeds
L
r
=
25 km. If the rupture, however, occurs along a smooth
(mature) fault, we observe a very strong increase of
L
r
with increas-
ing
D
p
≤
5.0 m (Fig. 7c, right-hand side). The PDFs hardly change
for
D
p
>
5.0 m, suggesting that the remaining rupture length is
largely independent from the present slip and every second rupture
will grow
L
r
>
65 km, that is,
M
>
7.0.
These findings imply that we might have a critical slip amplitude
(here
D
p
,
critical
∼
5.0 m) on a smooth (mature) fault, for which we
have an increased probability for the occurrence of a large earth-
quake for which EEW is needed. For the ruptures along heteroge-
neous (generic) faults, however, we cannot identify a critical slip
value (Fig. 7c, left-hand side). Note again, that
D
p
,
critical
as specified
here depends on the scaling of the 1-D slip models, and might not
be equal to the critical slip amplitude that we have to expect for
a real earthquake. A more accurate quantification of
D
p
,
critical
will
require further statistical analyses of earthquake data and modelling
of earthquake ruptures that is beyond the scope of this paper.
The statistical analyses of slip and rupture length in this paper
demonstrate that a large slip amplitude
D
p
does not necessarily mean
that we have a high probability for a large magnitude earthquake.
Or, a small-to -moderate slip amplitude does not necessarily mean
that the probability for a large event is low. We need to consider the
respective characteristics of the underlying fault. We conclude from
this finding, that an EEW system for large earthquakes requires
a mechanism for the rapid recognition of the causative fault and
consideration of its smoothness, that is, slip heterogeneity.
The majority of the algorithms, which currently are applied to
EEW, use a fixed time window of a couple of seconds of the seismic
C
2010 The Authors,
GJI
,
183,
1014–1030
Geophysical Journal International
C
2010 RAS
Early warning for large earthquakes
1023
Figure 10.
Probability of exceedance (POE) of different magnitude levels
M
for an ongoing rupture with a present magnitude of
M
p
5.5 and present slip
amplitude
D
p
. The probabilities are derived from the 1-D stochastic slip models for rough (left-hand side) and smooth ruptures (right-hand side), respectively.
The POE is controlled by the slip amplitudes,
D
p
, the smoothness of the causative fault, however, has an even larger impact. The probabilities for a large event
are much higher than in the G-R model (Fig. 1).
Table 1.
Probability of exceedance (POE) of different magnitude levels
M
for a rupture with a present
magnitude of
M
p
5.5. The POE is derived from 1-D stochastic slip models with two levels of rupture
roughness
μ
and slip amplitudes
D
p
. The probabilities are higher than predicted by the Gutenberg-and-Richter
(G-R) model.
Stochastic slip models with
μ
=
1.1/1.5 (rough/smooth rupture) and slip
D
p
Magnitude
G-R model
D
p
=
0.5 m
D
p
=
2.5 m
D
p
=
5m
D
p
=
10 m
D
p
=
15 m
M
≥
6.0
0.3
0.4/0.7
0.6/1.0
0.9/1.0
1.0/1.0
1.0/1.0
M
≥
6.5
0.1
0.2/0.5
0.3/0.8
0.5/0.9
0.8/1.0
0.9/1.0
M
≥
7.0
0.03
0.1/0.3
0.1/0.5
0.2/0.6
0.4/0.7
0.5/0.8
M
≥
7.5
<
0.01
<
0.1/0.1
<
0.1/0.2
0.1/0.2
0.1/0.3
0.1/0.3
P
wave for a rapid estimation of the earthquake magnitude, using,
for example, the average or the predominant period of shaking (e.g.
Nakamura 1988; Allen & Kanamori 2003; Kanamori 2005; Lewis
& Ben-Zion 2007; Wu
et al.
2007; Lewis & Ben-Zion 2008; B
̈
ose
et al.
2009a,b). This approach is heavily disputed, since the rupture
process of earthquakes with
M
>
6 takes much longer than 3 s and it
appears questionable, why the eventual rupture length and therewith
the magnitude of an earthquake should be pre-determined at this
early stage of the rupture process (e.g. Rydelek & Horiuchu 2006;
Rydelek
et al.
2007). Only if a large amount of the seismic energy of
the earthquake is released at its initial stage, that is, only if seismic
ruptures start in patches of large seismic slip (as proposed, e.g.
by Mai
et al.
2005 or Manighetti
et al
. 2005), a fast and robust
magnitude prediction might become feasible. If this is not the case,
the magnitude will likely remain underestimated, if we use a limited
time window of a few seconds only.
Our analyses might shed some light onto these discussions: it
seems appropriate to assign to a smooth (mature) fault, such as the
SAF, a higher probability that a rupture will continue and finally
evolve into a ‘Big One’, than to a heterogeneous (generic) fault. This
implies that the observation of a (apparent) moderate earthquake
(judging from the first few seconds of waveform data) on a smooth
(mature) fault should result in a stronger warning.
Once the probabilities of the total rupture length
L
and the even-
tual magnitude
M
of an ongoing earthquake are estimated (pro-
vided that
L
p
or
M
p
, respectively, are known), we can predict in a
probabilistic manner the level and distribution of ground shaking
amplitudes, such as PGV, along the evolving rupture using em-
pirical magnitude/distance relationships (see, e.g.
Next generation
attenuation
—special issue by
Earthquake Spectra
, 2008). The con-
sideration of rupture dimensions allows for much more realistic
assessments of the expected ground shaking than could be deliv-
ered by a point source approximation, in which the ground motion
amplitudes are assumed to decay with increasing epicentral rather
than rupture-to-site distance.
Peak ground displacement (PGD) [or inferred from this the re-
sponse spectral acceleration (SA) at long periods], on the other
hand, shows a good correlation with the static slip on the fault (e.g.
Aagaard
et al.
2001). This means that if we were able to predict the
future slip of the evolving rupture, we could use this information
for a significant improvement of the estimates of parameters PGD
and SA for EEW.
Demonstration for an
M
7.8 scenario earthquake
along the southern SAF
We are going to demonstrate the probabilistic prediction of seismic
ground motion and future slip for an
M
7.8 scenario earthquake
along the southern SAF in California (Fig. 11a, left-hand side). The
SAF is a right-lateral strike-slip fault with a total length of about
1200 km. The average time interval between significant ruptures
along the SAF varies from about 25 yr in the Cholame Valley to
C
2010 The Authors,
GJI
,
183,
1014–1030
Geophysical Journal International
C
2010 RAS