Geophysical Journal International
Geophys. J. Int.
(2012)
191,
803–812
doi: 10.1111/j.1365-246X.2012.05657.x
GJI Seismology
Real-time Finite Fault Rupture Detector (FinDer)
for large earthquakes
Maren B
̈
ose, Thomas H. Heaton and Egill Hauksson
Seismological Laboratory, California Institute of Technology,
1200
E. California Blvd., Mail Code
252
–
21
, Pasadena,
CA 91125
,USA.
E-mail: mboese@caltech.edu
Accepted 2012 August 20. Received 2012 June 29; in original form 2012 January 10
SUMMARY
To provide rapid estimates of fault rupture extent during large earthquakes, we have developed
the
Fin
ite Fault Rupture
De
tecto
r
algorithm, ‘FinDer’. FinDer uses image recognition tech-
niques to detect automatically surface-projected fault ruptures in real-time (assuming a line
source) by estimating their current centroid position, length
L
, and strike
θ
. The approach is
based on a rapid high-frequency near/far-source classification of ground motion amplitudes in
a dense seismic network (station spacing
<
50 km), and comparison with a set of pre-calculated
templates using ‘Matching by Correlation’. To increase computational efficiency, we perform
the correlation in the wavenumber domain. FinDer keeps track of the current dimensions of
a rupture in progress. Errors in
L
are typically on the same order as station spacing in the
network. The continuously updated estimates of source geometries as provided by FinDer
make predicted shaking intensities more accurate and thus more useful for earthquake early
warning, ShakeMaps, and related products. The applicability of the algorithm is demonstrated
for several recorded and simulated earthquakes with different focal mechanisms, including the
2009
M
w
6.3 L’Aquila (Italy), the 1999
M
w
7.6 ChiChi (Taiwan) and the
M
w
7.8 ShakeOut
scenario earthquake on the southern San Andreas Fault (California).
Key words:
Image processing; Spatial analysis; Earthquake ground motions; Seismic mon-
itoring and test-ban treaty verification; Earthquake interaction, forecasting and prediction;
Early warning.
INTRODUCTION
Using current data acquisition, telemetry and processing technolo-
gies, seismologists are able to determine within seconds the mag-
nitude and location of a local earthquake using dense networks of
seismic sensors with real-time communication to a central process-
ing site. However, the current processing algorithms are insufficient
in the case of large earthquakes (
M
w
>
6); large earthquakes rupture
for long distances along seismic faults (
M
w
7
∼
60 km,
M
w
7.5
∼
130 km, and
M
w
8
∼
300 km; Wells & Coppersmith 1994); it is the
distance to the closest fault rupture that causes the peak strength of
shaking at a given site. To estimate the distribution of shaking and
potential damage, for example, for earthquake early warning (EEW)
in the seconds before seismic waves hit or for the coordination of
emergency responders afterwards, knowledge on the approximate
location, length and orientation of the earthquake rupture is crucial.
Estimating fault rupture extent in real-time though is difficult; the
temporal and spatial evolution of extended fault ruptures can be
highly complex.
The majority of EEW systems estimate the expected shaking
intensity at a given user site from the magnitude and epicentral
distance
R
epi
of the earthquake applying empirical ground motion
prediction equations (e.g. Allen
et al.
2009; B
̈
ose
et al.
2012a). For
larger earthquakes (
M
w
>
6), however, it would be more appropriate
to use rupture-to-site distances, such as Joyner–Boore distance
R
jb
,
which required knowledge on source finiteness (B
̈
ose & Heaton
2010). Neglecting the rupture extent can underestimate shaking
intensities and result in warnings not being issued because shaking
predictions do not meet predefined thresholds.
A prominent example is from the recent 2011
M
w
9 Tohoku-
Oki earthquake in Japan; while ground motion predictions by the
Japanese Meteorological Agency (JMA) were fairly accurate in
the Sendai area, which is close to the point of rupture nucleation,
a warning was not transmitted to people and users in the Kanto
region around Tokyo. Even though Kanto shaking was strong, it
was so far from the epicentre (
R
epi
>
350 km) that it was seriously
under predicted by the JMA EEW system (Hoshiba & Iwakiri 2011;
Sagiya
et al.
2011). The recognition of fault rupture extent would
likely have significantly improved the system performance during
the Tohoku-Oki earthquake (Kurahashi & Irikura 2011).
Information on earthquake source finiteness is also important for
the computation of rapid response ShakeMaps, which are maps of
spatially interpolated recorded ground motions. These maps pro-
vide critical information to search and rescue teams and other
C
2012 The Authors
803
Geophysical Journal International
C
2012 RAS
804
M. B
̈
ose, T.H. Heaton and E. Hauksson
stakeholders in the aftermath of strong earthquakes (Wald
et al.
1999). Generally, fault finiteness is not explicitly considered in
ShakeMaps within the first hour(s) following the earthquake. In
the case of the 2011 Tohoku-Oki earthquake, for instance, a first
ShakeMap considering source finiteness was available 2–3 hr af-
ter the event (Hayes
et al.
2011). Currently, fault dimensions in
ShakeMaps are estimated from near-fault strong motion data, sur-
face offset observations, geodetic displacements, regional and lo-
cal waveforms and teleseismic data, which need to be processed
and interpreted (Wald
et al.
2005). Recently, Convertito
et al.
(2012) proposed a Bayesian approach to estimate fault rupture ex-
tent within minutes using the residuals of observed and predicted
ground-motion quantities.
To provide estimates of fault rupture extent during large earth-
quakes while still in progress, we have developed the
Fin
ite Fault
Rupture
De
tecto
r
(FinDer) algorithm. FinDer uses image recogni-
tion techniques to detect and map finite fault ruptures automatically
and in real-time based on a rapid near/far-source classification of
seismic observations recorded by a dense seismic network (station
spacing
<
50 km), and comparison of these observations with a set
of pre-calculated templates. The continuously updated estimates of
source dimensions as provided by FinDer make predicted ground
motions more accurate (since they can be predicted from
R
jb
rather
than
R
epi
) and thus more useful for EEW, ShakeMaps and related
products.
FINITE FAULT RUPTURE DETECTOR
‘FinDer’
The Finite Fault Rupture Detector algorithm FinDer detects and
maps finite source earthquake ruptures (
M
w
>
6) in real-time. Given
data available at time
t
i
, the algorithm estimates the current centroid
position, rupture length
L
and strike
θ
of the assumed line source
(without weighting for slip). The algorithm accomplishes this by
comparing recorded ground motion amplitudes with a set of pre-
calculated templates. The processing flow is shown in Fig. 1.
FinDer presumes the existence of a dense seismic network (sta-
tion spacing
<
50 km) with real-time data transmission to a central
processing site. Let us assume that at a given time
t
i
a set of ground
motion observations at all or at a subset of seismic sensors is avail-
able. Whenever an earthquake is detected [for example, using the
Virtual Seismologist (Cua & Heaton 2007), PreSEIS (B
̈
ose
et al.
2008; B
̈
ose
et al.
2012b) or any other algorithm for EEW (Allen
et al.
2009)], ground motion amplitudes are soil-to-rock corrected,
Figure 1.
Processing steps of the Finite Fault Rupture Detector (FinDer) algorithm to determine the current centroid position, length
L
,andstrike
θ
of an
ongoing earthquake rupture. Step 1: gathering of maximum ground motion amplitudes in a seismic network at time
t
i
. Step 2: empirical site-corrections,
map projection, and spatial interpolation of amplitudes using zero-amplitude boundary condition. Step 3: near/far source classification, for exam
ple, based on
acceleration thresholds. Steps 4–6: comparison of binary map
f
(
x
,
y
) obtained from classification with a set of pre-calculated templates for various
L
and
θ
using ‘Matching by Correlation’. The best match is achieved at the point of maximum correlation (Step 6). For increased speed, correlation is perform
ed in the
wavenumber domain. Step 7: the template (and thus
L
and
θ
) with the best fit is obtained by minimizing the misfit (SSE error) between template and ground
motion observations using a Direct Pattern Search algorithm. Steps 1–7 are repeated regularly, for example, every
t
=
1to2s.
C
2012 The Authors,
GJI
,
191,
803–812
Geophysical Journal International
C
2012 RAS
Finite Fault Rupture Detector (FinDer)
805
Figure 2.
Saturation of peak ground acceleration (PGA) as a function of
moment magnitude
M
w
and Joyner–Boore distance
R
jb
using relations by
Cua (2005). A simple acceleration threshold is effective for near/far-source
classification, because large PGA values are typically expected for stations
close to the rupture, and these large PGA values occur for all events
M
w
>
6.
In this study, we apply three acceleration thresholds at 95, 70 and 55 cm s
–
2
,
to identify near-source stations for
M
w
>
6at
R
jb
≤
15,
≤
20 and
≤
25 km,
respectively.
projected onto a map, and interpolated spatially onto a fine grid
assuming a zero-amplitude boundary condition (Fig. 1, Steps 1
and 2). Shortly after rupture nucleation the seismic
P
and
S
waves
will have reached only a few sensors in the network. The spatial
interpolation of maximum amplitudes is carried out over all sta-
tions, independent from the type of observation (
P
wave,
S
wave,
surface wave, background noise); that is, independent from wave
propagation.
In the next step, FinDer assigns a unit less value of ‘
+
1’ to all
sites in the interpolated map that appear to be close to the rupture,
and a value of ‘0’ to those further away. This near/far-source classi-
fication gives a binary image
f
(
x
,
y
) with Cartesian coordinates
x
and
y
(Fig. 1, Step 3). In this study, we use simple high-frequency (ac-
celeration) thresholds for near/far-source classification, motivated
by the observation that high acceleration values are typically ob-
served very close to the rupture only (Fig. 2). We apply acceleration
thresholds at 95, 70 and 55 cm s
–
2
, to identify near-source stations
for earthquakes with
M
w
≥
6 at Joyner–Boore distances of
R
jb
≤
15,
≤
20 km and
≤
25 km, respectively (Fig. 2). Principally, however,
we could also apply more sophisticated algorithms for classifica-
tion that make use, for example, of additional information from
mid- and low-frequency observations (Yamada
et al.
2007; B
̈
ose
et al.
2012b).
To estimate the current centroid position, length
L
,andstrike
θ
of the underlying rupture, FinDer compares the binary
image f
(
x
,
y
)
obtained from near/far-source classification with a set of around
12,500 pre-calculated
templates g
(
x
,
y
|
L,
θ
) for various
L
and
θ
(here:
5km
≤
L
≤
350 km with increments of
L
=
5 km, and 0
◦
≤
θ<
180
◦
with
θ
=
1
◦
). We assume that the spatial coverage of each
template is much smaller than that of the interpolated map. For a
systematic comparison of the impact of different near/far-source
distances in this paper, we generated three sets of templates. Sites
close to the rupture (
R
jb
≤
15,
≤
20 and
≤
25 km, respectively) were
assigned a unit less value of ‘
+
1’, and a value of ‘0’ to those further
away. To determine the optimum solutions of
L
and
θ
at time
t
i
,we
need to find (1) for a given
L
–
θ
pair, the place in
f
(
x
,
y
) that shows
the best match with template
g
(
x
,
y
|
L,
θ
), and (2) the template with
the best overall fit, that is smallest error.
To solve requirements (1) and (2) in real-time, we use the tem-
plates g(
x
,
y
|
L,
θ
) as spatial filters and compute the sum of products
for each location of
g
in
f
. The best match(es) is (are) the loca-
tion(s) of the maximum value(s) in the resulting correlation image.
This method is known as ‘Matching by Correlation’ (e.g. Gonzales
et al.
2004). However, since this approach is generally computation-
ally intensive and therefore not well suited for real-time procedures
such as EEW, we implement correlation in the wavenumber domain,
using the correlation theorem:
f
(
x
,
y
)
g
(
x
,
y
|
L
,θ
)
⇔
̃
f
(
k
x
,
k
y
)
̃
g
∗
(
k
x
,
k
y
|
L
,θ
)
.
(1)
Eq. (1) relates the spatial correlation to the product of the Fourier
transforms
̃
f
(
k
x
,
k
y
)and
̃
g
∗
(
k
x
,
k
y
|
L
,θ.
) in the wavenumber do-
main (
k
x
,
k
y
), where ‘
’ denotes correlation and ‘
∗
’ the complex
conjugate. In other words, spatial correlation is obtained as the
inverse Fourier transform of the product of the transform of one
function times the conjugate of the transform of the other (e.g.
Gonzales
et al.
2004). This approach is very efficient and can be
parallelized for additional speed if needed.
The maximum of the correlation image
̃
f
(
k
x
,
k
y
)
̃
g
∗
(
k
x
,
k
y
|
L
,θ
)
gives the location of the best match of
subimage f’
(
x
,
y
) and template
g(
x
,
y
|
L,
θ
)foragiven
L–
θ
pair (Fig. 1, Steps 4–6). If we suspect that
multiple earthquakes occur simultaneously, we can analyse further
(local) maxima in the correlation image.
To determine the optimum solutions of
L
and
θ
,
we need to find
the template with the smallest misfit, that is, the minimum sum of
squared errors (SSE):
SSE
=
∑
x
∑
y
[
f
(
x
,
y
)
−
g
(
x
,
y
|
L
,θ
)
]
.
(2)
We use a Direct Pattern Search optimization algorithm (Math-
works Matlab Global Optimization Toolbox), which usually finds
a minimum of eq. (2) after less than 15 iterations. Thus we can
quickly estimate the current location and dimensions of the on-
going rupture, and determine
L
and
θ
in real-time (Fig. 1, Step
7). Once the current rupture length
L
is determined
,
we can also
estimate magnitude
M
w
from empirical
L–M
w
relations, such as
proposed by Wells & Coppersmith (1994). Note, however, that
these relations depend on the type of focal mechanism and on
whether the rupture is along the surface or subsurface, involv-
ing additional uncertainties of
∼±
0.1–0.2
M
w
units in the
L
–
M
w
relationship.
Steps 1–7 are repeated every
t
seconds using the currently avail-
able ground motion amplitudes. The total processing time depends
on the desired resolution and spatial coverage of the interpolated
maps and templates. In this study, we used a resolution of 0.025
◦
for
L
≤
50 km and a resolution of 0.05
◦
for
L
>
50 km. For the
examples shown in this paper, the computation of
L
and
θ
at each
time
t
i
takes on the order of around 1 s on a modern desktop com-
puter, that is, we can update the estimates every
t
=
1–2 s. A
higher resolution would increase the computational effort without
gaining much additional information about the rupture of a large
earthquake.
DATA
In this study, we apply FinDer to several recorded and simu-
lated earthquakes with different focal mechanisms, including the
2009
M
w
6.3 L’Aquila (Italy) normal fault earthquake (e.g. Ameri
et al.
2012; Gallovi
ˇ
c&Zahradn
́
ık 2012), the 1999
M
w
7.6 ChiChi
C
2012 The Authors,
GJI
,
191,
803–812
Geophysical Journal International
C
2012 RAS
806
M. B
̈
ose, T.H. Heaton and E. Hauksson
Table 1.
Observed and estimated rupture length
L
obs
and
L
est
, and azimuth
θ
obs
and
θ
est
at six time steps after
origin. Errors in
L
are on the same order as station spacing. Map resolution is 0.05
◦
(
∼
5 km), that is
L
est
can change
in these examples in increments of
L
=
5 km only.
L
est
,
15
and
θ
est
,
15
were estimated using a near/far-source
classification threshold of 95 cm s
–
2
(
R
jb
≤
15 km);
L
est
,
20
and
θ
est
,
20
refer to a near/far-source classification
thresholdof70cms
–
2
(
R
jb
≤
20 km), and
L
est
,
25
and
θ
est
,
25
refer to 55 cm s
–
2
(
R
jb
≤
25 km). See Figs 4 and 5
for plotted data. Note that the rupture duration of the
M
7.8 ShakeOut scenario earthquake exceeds 60 s.
Event
10 s
20 s
30 s
40 s
50 s
60 s
Final
M
6.3 L’Aquila (RAN)
L
obs
20 km
20 km
20 km
20 km
20 km
20 km
20 km
L
est
,
15
5km
10km
15km
20km
20km
20km
20km
L
est
,
20
5km
10km
10km
10km
10km
10km
10km
L
est
,
25
5km
15km
20km
30km
30km
30km
30km
θ
obs
140
◦
140
◦
140
◦
140
◦
140
◦
140
◦
140
◦
θ
est
,
15
45
◦
66
◦
84
◦
91
◦
110
◦
110
◦
110
◦
θ
est
,
20
135
◦
133
◦
136
◦
136
◦
139
◦
139
◦
139
◦
θ
est
,
25
133
◦
127
◦
120
◦
117
◦
110
◦
110
◦
110
◦
M
7.6 ChiChi (TSMIP)
L
obs
40 km
80 km
100 km
100 km
100 km
100 km
100 km
L
est
,
15
5 km
45 km
80 km
110 km
120 km
120 km
120 km
L
est
,
20
10 km
45 km
80 km
110 km
125 km
135 km
135 km
L
est
,
25
5 km
35 km
75 km
105 km
125 km
130 km
130 km
θ
obs
23
◦
23
◦
23
◦
23
◦
23
◦
23
◦
23
◦
θ
est
,
15
45
◦
37
◦
37
◦
37
◦
34
◦
36
◦
36
◦
θ
est
,
20
132
◦
73
◦
51
◦
42
◦
23
◦
18
◦
18
◦
θ
est
,
25
43
◦
20
◦
28
◦
33
◦
37
◦
46
◦
46
◦
M
7.8 ShakeOut (CISN)
L
obs
50 km
80 km
110 km
135 km
165 km
195 km
300 km
L
est
,
15
5 km
40 km
70 km
110 km
145 km
225 km
335 km
L
est
,
20
5 km
45 km
85 km
115 km
160 km
225 km
330 km
L
est
,
25
10 km
50 km
90 km
125 km
165 km
240 km
350 km
θ
obs
121
◦
121
◦
121
◦
121
◦
121
◦
121
◦
121
◦
θ
est
,
15
101
◦
120
◦
120
◦
120
◦
120
◦
120
◦
120
◦
θ
est
,
20
146
◦
125
◦
123
◦
122
◦
119
◦
119
◦
119
◦
θ
est
,
25
139
◦
127
◦
124
◦
122
◦
119
◦
119
◦
119
◦
(Taiwan) thrust fault earthquake (e.g. Ji
et al.
2001) and the 2008
M
w
7.8 ShakeOut scenario (southern California) strike-slip fault
earthquake (Jones
et al.
2008).
On 2009 April 6, a
M
w
6.3 earthquake hit the city of L’Aquila
(Central Italy), located at only a few kilometres northeast of the
epicentre, and causing more than 300 casualties and large dam-
age to the city and neighboured villages (Ameri
et al.
2012). The
earthquake occurred along a NW–SE trending normal fault, with an
estimated rupture length of
L
obs
≈
20 km and strike of
θ
obs
≈
140
◦
(Ameri
et al.
2012). We downloaded the corrected strong-motion
waveforms recorded by the Italian Strong motion Network (
Rete
Accelerometrica Nazionale
, RAN) from the Italian strong-motion
database (ITACA, http://itaca.mi.ingv.it); we then determined the
maximum absolute ground motion amplitudes (largest value of all
three components) of the corrected records in time windows of 1 s
on each record to obtain acceleration envelopes, such as described
in Cua & Heaton (2007). About 14 of the waveforms were recorded
within 50 km distance from the earthquake epicentre (Ameri
et al.
2012).
The 1999
M
w
7.6 ChiChi earthquake (Taiwan) was one of the
largest inland events in the 20th century and produced a rela-
tively dense set of strong motion recordings. It produced an almost
100-km-long surface rupture along the primarily NS striking
Chelungpu fault with offsets of up to 8 m near the northern end of
the fault (Ji
et al.
2001). The Harvard CMT solution shows a shallow
(27
◦
) dipping thrust movement towards east. Due to this small dip
angle we do not expect the surface-projected rupture (as predicted
by FinDer) to agree with the surface trace of the Chelungpu fault,
but to be shifted to the east. From a fault model proposed by Ji
et al.
(2001) we estimate the rupture length, azimuth, and velocity
of the ChiChi earthquake as
L
obs
≈
100 km,
θ
obs
≈
23
◦
,and
v
r
=
2.0 km s
–
1
, respectively (Table 1).
Seismic waveforms of the ChiChi earthquake were recorded at
441 free-field strong-motion stations of the Taiwan Strong-Motion
Instrumentation Program (TSMIP) operated by the Central Weather
Bureau (CWB) with an average station spacing of
∼
5 km. We
downloaded the acceleration records from the COSMOS Virtual
Datacenter (http://db.cosmos-eq.org), and determined the maxi-
mum absolute ground motion amplitudes of the corrected records
in time windows of 1 s of each record to obtain acceleration en-
velopes (Fig. 3a). We apply Vs
30
values from the National Cen-
ter for Research on Earthquake Engineering (NCREE) and CWB
(http://geo.ncree.org.tw) for empirical soil-to-rock corrections (ex-
planation follows).
The ShakeOut scenario earthquake is a best-estimate scenario
rupture of a
M
w
7.8 earthquake on the southernmost portion of
the San Andreas Fault (Jones
et al.
2008). The rupture of this right-
lateral strike slip event starts at the Salton Sea in the south and travels
over a distance of
L
obs
=
300 km in a northwesterly direction towards
Lake Hughes. The ShakeOut scenario earthquake was used to study
C
2012 The Authors,
GJI
,
191,
803–812
Geophysical Journal International
C
2012 RAS
Finite Fault Rupture Detector (FinDer)
807
Figure 3.
Examples for (a) observed and (b) simulated acceleration envelope functions for the
M
w
7.6 ChiChi and the
M
w
7.8 ShakeOut scenario earthquake
rupture (black lines). The simulations in (b) were obtained by summing the envelopes of multiple smaller subsource events of
M
w
6 (grey stars). We use
high-frequency thresholds at 95, 70 and 55 cm s
–
2
, to identify near-source stations for earthquakes with
M
w
>
6 at Joyner–Boore distances of
R
jb
≤
15,
≤
20
and
≤
25 km, respectively (here we show classification results for 55 cm s
–
2
only). The envelopes are used to determine if and when these thresholds at a given
station are reached.
the potential impact of a major earthquake on southern California
and also formed the bases for the 2008 November 13, ShakeOut
Earthquake Drill with over five million participants (Jones
et al.
2008).
Graves
et al.
(2010) simulated seismic ground-motions for the
ShakeOut scenario earthquake in the low- to mid-frequency band
using a finite difference method. However, these waveforms cannot
be used for the testing of FinDer, since our near/far-source clas-
sification is based on high-frequency motions. We therefore apply
a method proposed by Yamada & Heaton (2008) to simulate the
envelope time functions of high-frequency ground motions during
large earthquakes.
Simulated envelope time series
To calculate the envelope time functions of ground motions for the
ShakeOut scenario earthquake, we apply a simple 2D source-model
C
2012 The Authors,
GJI
,
191,
803–812
Geophysical Journal International
C
2012 RAS
808
M. B
̈
ose, T.H. Heaton and E. Hauksson
approach (Yamada & Heaton 2008) that extends earlier ground-
motion models by Cua (2005) and Cua & Heaton (2007). In this
approach the fault surface is divided into multiple (not overlapping)
subfaults, each represented by a single point (sub-) source that
radiates
P
and
S
waves when the rupture front arrives. The envelope
of the
i
th subsource is described by
E
i
(
t
∣
∣
M
,
R
epi
)
=
E
[
t
−
t
i
∣
∣
M
i
,
R
epi
,
i
]
,
(3)
where
M
i
is the magnitude of the subsource
i
,
R
epi
,
i
is the (epicentral)
subsource-to-site distance, and
t
i
is the time delay due to the rupture
propagation. Eq. (3) can be solved analytically for point source
earthquakes (
M
≤
6) from ground-motion models proposed by Cua
(2005) and Cua & Heaton (2007).
For larger finite-source earthquakes (
M
>
6), the total envelope
of ground motion
E
total
at a given site is modelled as the combination
of the responses of each subsource
i
(Fig. 3b). For high-frequency
motions with roughly random phase, the square root of the sum of
the squares of the envelope amplitudes from each subsource gives
a good approximation of
E
total
:
E
total
(
t
∣
∣
M
,
R
jb
)
=
√
∑
N
i
=
0
E
2
i
(
t
∣
∣
M
,
R
epi
)
,
(4)
where
N
is the number of the point sources. This approach is only
suited for the modelling of high-frequency ground motions that
seem to be fairly insensitive to the effects of source radiation and
directivity (Yamada & Heaton 2008).
We assume that the dimensions of all subsources are uniform
and that each subsource magnitude is
M
i
=
6, corresponding to a
subfault length of
l
≈
10 km (Wells & Coppersmith 1994). The
number of subsources is determined from
N
≈
(
L
/
l
). We use a 1-D
seismic velocity model for southern California to determine the
P
and
S
wave arrivals for the ShakeOut scenario (Hutton
et al.
2010).
Fig. 3(b) shows three example simulations.
Cua (2005) and Cua & Heaton (2007) determined ground motion
models for both rock and soil conditions. They classified sites with
V
s30
≤
434 m s
–
1
as ‘soil’ and sites with
V
s30
>
434 m s
–
1
as
‘rock’, where
V
s30
is the average shear wave velocity taken over the
uppermost 30 m. We simulate a set of envelope functions for both
rock and soil condition for
M
6 and distances of up to 100 km, and
obtain a mean amplification factor of peak amplitudes in the order
of 1.35. We take this empirical value for soil-to-rock corrections of
simulated and observed ground motion envelopes in this study.
We simulate the rupture of the ShakeOut scenario earthquake
propagating with a rupture speed of
v
r
= 2.9 km s
–
1
, using the cur-
rent distribution of seismic real-time stations in the southern Cal-
ifornia Integrated Seismic Network (CISN/SCSN; www.cisn.org
and www.scsn.org) with an average station spacing of
∼
20 km.
The rupture of the
M
w
7.8 ShakeOut scenario earthquake extends
across the Big Bend of the San Andreas Fault. We approximate this
trace by a linear rupture with
θ
obs
=
121
◦
, which requires shifting
the epicentre by 0.15
◦
(
∼
17 km) to the east (Table 1). We apply
V
s30
values determined by Wills
et al.
(2000) for the simulation of
envelope functions at rock and soil stations.
Alike for the recorded L’Aquila and ChiChi earthquakes, we
apply thresholds of 95, 70 and 55 cm s
–
2
to the simulated ground-
motion envelopes to estimate, whether a site is close to the rupture
(
≤
15,
≤
20 or
≤
25 km, respectively) or further away. The simula-
tions help us to identify if and when these thresholds at a given
station are exceeded.
RESULTS
We replay the records of each earthquake as simulated real-time
data streams and emulate the temporal evolution of predictions by
FinDer at regular time steps
t
i
after origin time (O.T.). We assume
that the transmission delays of seismic waveform data from the
stations to the central processing facility are insignificantly small
and negligible.
Table 1 and Fig. 4 compare the observed and estimated rupture
lengths and azimuths,
L
obs
,
L
est
,
θ
obs
and
θ
est
, for each of the three
earthquakes. Values for
L
obs
and
θ
obs
are taken from the literature
and were discussed earlier;
L
est
,
15
and
θ
est
,
15
were estimated using a
near/far-source classification threshold of 95 cm s
–
2
(
R
jb
≤
15 km);
L
est
,
20
and
θ
est
,
20
refer to a near/far-source classification threshold
of 70 cm s
–
2
(
R
jb
≤
20 km), and
L
est
,
25
and
θ
est
,
25
refer to 55 cm s
–
2
(
R
jb
≤
25 km). Fig. 5 shows screenshots of the estimated near/far
source station classifications (using an acceleration threshold of
70 cm s
–
2
) and estimated rupture parameters in map view at four
time steps after origin.
For all events we observe a good convergence towards the ob-
served values
L
obs
and
θ
obs
with increasing time. The smaller the
station spacing, the quicker the detection of a large earthquake, and
the more accurate are the estimated parameters (Table 1; Figs 4 and
5). The map resolution in these examples is 0.05
◦
(
∼
5 km), that
is
L
est
increases in increments of 5 km only. Note that the earth-
quake ruptures are usually still in progress when FinDer provides
the first estimates of their dimensions and azimuth. Usually, the
predictions of
L
are smaller or equal to the current rupture length;
that is, FinDer keeps track of the current rupture parameters without
predicting their future evolution (Table 1; Figs 4 and 5). Generally,
errors in
L
are on the same order as the station spacing. There are no
systematic differences in the performance depending on the three
near/far-source classification thresholds.
The rupture length of the ChiChi earthquake appears to be over-
estimated by
∼
30 per cent around 1 min after rupture nucleation
(Table 1, Fig. 5b). Shin & Teng (2001), however, suggested the oc-
currence of a secondary rupture along the Shihtan fault shortly after
the ChiChi main shock. This rupture is consistent with an extended
zone of near-source classified stations, north of the main rupture, as
was also observed by Yamada
et al.
(2007). Both ruptures occur very
close to each other in space and time and cannot be distinguished
by FinDer. Also note that during the ChiChi earthquake some sta-
tions along the Taiwanese east coast were incorrectly classified as
near-source sites (Fig. 5b); the predictions of
L
and
θ
,however,are
robust and not affected by these misclassifications.
Estimated shaking intensity
Most EEW systems, including the current version of the CISN
ShakeAlert Earthquake Early Warning Demonstration System for
California (B
̈
ose
et al.
2012a), estimate the shaking intensity at a
given user site from predicted magnitudes and epicentral distances
R
epi
using empirical ground motion prediction equations. More ac-
curate predictions are obtained, if
R
epi
is replaced by Joyner–Boore
distance
R
jb
; that is, the closest distance of the user site to the
surface-projected earthquake rupture (B
̈
ose & Heaton 2010). This
procedure, however, requires the knowledge or estimation of rupture
dimensions, such as those provided by FinDer.
In the following, we will analyse the estimated and observed
Modified Mercalli Intensity (MMI) intensities for the
M
w
7.8 Shake-
Out scenario earthquake in three cities in southern California: Palm
Springs (
R
epi
=
80 km), Riverside (
R
epi
=
160 km) and Los Angeles
C
2012 The Authors,
GJI
,
191,
803–812
Geophysical Journal International
C
2012 RAS
Finite Fault Rupture Detector (FinDer)
809
Figure 4.
Comparison of observed and estimated rupture length
L
and strike
θ
for the
M
w
6.3 L’Aquila,
M
w
7.6 ChiChi and
M
w
7.8 ShakeOut scenario
earthquakes within the first 2 min after origin. The dashed lines show the current rupture length assuming constant rupture velocities of
v
r
=
2.0 km s
–
1
for
the ChiChi earthquake, and
v
r
=
2.9 km s
–
1
for the two other events. FinDer keeps track of the current rupture parameters; estimates in
L
are usually below
or equal to the current rupture dimensions. The plots show the results for three templates with near/far-source classification thresholds at 95 cm s
–
2
(
R
jb
≤
15 km), 70 cm s
–
2
(
R
jb
≤
20 km) and 55 cm s
–
2
(
R
jb
≤
25 km).
Figure 5.
Results of near/far-source classification and predicted (surface-projected) ruptures for the (a)
M
w
6.3 L’Aquila (Italy), (b)
M
w
7.6 ChiChi (Taiwan),
and (c)
M
w
7.8 ShakeOut scenario (southern California) earthquakes at four time steps after rupture nucleation. Black squares mark near-source, grey triangl
es
far-source classified stations determined from high-frequency thresholds (here: 70 cm s
–
2
). Circles show the
P
and
S
wave fronts; grey lines show evolution of
L
obs
. The grey rectangles in (b) show a three-plane fault geometry with shallow dip towards east determined from static inversion (Ji
et al.
2001). The estimated
ruptures are robust and not strongly affected by misclassified stations in (b). The final fault rupture extent of the ChiChi earthquake is overestimate
dby
∼
30 per cent, caused by a secondary rupture along the Shihtan fault in the north shortly after the nucleation of the main shock.
C
2012 The Authors,
GJI
,
191,
803–812
Geophysical Journal International
C
2012 RAS