of 13
Real-Time Estimation of Fault Rupture Extent
Using Envelopes of Acceleration
by Masumi Yamada and Thomas Heaton
Abstract
We present a new strategy to estimate the geometry of a rupture on a
finite fault in real time for earthquake early warning. We extend the work of Cua
and Heaton who developed the virtual seismologist (
VS
) method (Cua, 2005), which
is a Bayesian approach to seismic early warning using envelope attenuation relation-
ships. This article extends the
VS
method to large earthquakes where fault finiteness is
important. We propose a new model to simulate high-frequency motions from earth-
quakes with large rupture dimension: the envelope of high-frequency ground motion
from a large earthquake can be expressed as a root-mean-squared combination of en-
velope functions from smaller earthquakes. We use simulated envelopes of ground
acceleration to estimate the direction and length of a rupture in real time. Using
the 1999 Chi-Chi earthquake dataset, we have run simulations with different para-
meters to discover which parameters best describe the rupture geometry as a function
of time. We parameterize the fault geometry with an epicenter, a fault strike, and two
along-strike rupture lengths. The simulation results show that the azimuthal angle of
the fault line converges to the minimum uniquely, and the estimation agrees with the
actual Chi-Chi earthquake fault geometry quite well. The rupture direction can be
estimated at 10 s after the event onset, and the final solution is achieved after
20 s. While this methodology seems quite promising for warning systems, it only
works well when there is an adequate distribution of near-source stations.
Introduction
Recently, with advances in data analysis and increased
awareness of the seismic hazard, the topic of earthquake
early warning has attracted more research attention and var-
ious early warning methods have been proposed from seis-
mologists and engineers (Nakamura, 1988; Allen and
Kanamori, 2003; Odaka
et al.
, 2003; Wu and Kanamori,
2005). Currently, the most ambitious system is the earth-
quake early warning system provided by the Japan Meteor-
ological Agency, which was released to the public in October
2007. The news of the system was broadcast widely, and it
attracted considerable public attention in Japan. The goal of
seismic early warning is to initiate optimal mitigating actions
based on the arrival time and amplitude of seismic waves
predicted at a given location. To achieve this, an earthquake
early warning system must collect and quickly analyze seis-
mic data in a manner that can be used to predict future shak-
ing. In principle, this could be achieved by using the present
value of an approximately known wave field as a boundary
condition to predict future wave fields using Navier
s equa-
tion (Baker
et al.
, 2005). However, from a practical view-
point, there are advantages to data analysis schemes that
involve characterization of the earthquake source. Predic-
tions of future shaking can be achieved by utilizing the ex-
tensive existing work on predicting ground shaking from
seismic sources. Ideally, an early warning system would pro-
vide the best estimate of slip in time and space that can be
deduced from seismic data available at any given instant in
time. We introduce a two-step strategy to accomplish this:
(1) We determine the spatial and temporal extent of an on-
going rupture by analyzing waveform envelopes of high-
frequency shaking, and (2) we determine approximate slip
from simple projections of long-period shaking onto the ap-
proximately known location of the rupture. It seems likely
that true ground displacement may be very helpful in the
real-time determination of fault slip. Because it is difficult
to obtain displacement from inertial seismometers (Clinton,
2004), the future incorporation of real-time high-sample-rate
Global Positioning System (
GPS
) into early warning systems
may be quite important. A more complete description of our
strategy to determine distributions in real time is given by
Yamada
et al.
(2007).
Cua and Heaton developed the virtual seismologist (
VS
)
method (Cua, 2005). It is a Bayesian approach to seismic
early warning designed for modern seismic networks and
proposed for small to moderate earthquakes with ruptures
that can be approximately modeled as a point source. The
607
Bulletin of the Seismological Society of America, Vol. 98, No. 2, pp. 607
619, April 2008, doi: 10.1785/0120060218
VS
algorithm uses an envelope attenuation relationship and
the predominant frequency content from the first few seconds
after the
P
-wave arrival. The advantage of the
VS
method is
its capacity to assimilate different types of information that
may be useful to find quick and reliable estimates of mag-
nitude and location (Cua, 2005). It gives the best estimate of
an earthquake property in terms of a probability density func-
tion. The Bayesian approach is a scheme to emulate human
capabilities to judge complex information by modeling un-
certainty in a probabilistic way.
Our goal is to extend the
VS
method to large earthquakes
where fault finiteness is important. Although there are some
researchers trying to characterize fault finiteness in real time
(Nakamura, 2001), most of earthquake early warning sys-
tems focus on estimating epicenters and magnitudes of earth-
quakes, not the fault geometry (Allen and Kanamori, 2003).
However, for large earthquakes, rupture length can be on the
order of tens to hundreds of kilometers, and the inhomoge-
neous slip distribution significantly affects the ground-
motion amplitude at a site. For example, the fault rupture
in the 1999 Chi-Chi earthquake was longer than 80 km,
and the largest slip was recorded at the northern end of
the fault. It would be difficult, if not impossible, to predict
such large shaking at large distances from the epicenter when
using a scheme that only characterizes the earthquake as a
point source.
In this article, we introduce a methodology to determine
the direction and length of the rupture from high-frequency
ground motions. We construct a model to describe high-
frequency motions from earthquakes with large fault dimen-
sion. By applying this model to high-frequency ground
motions in real time, it becomes possible to estimate the geo-
metry of an ongoing earthquake.
Data
The data for this analysis is the strong-motion dataset
from the 21 September 1999 Chi-Chi Earthquake that oc-
curred in central Taiwan (Lee
et al.
, 2001). The epicenter
was located at 120.82° N, 23.85° E, with a focal depth of
8 km, according to the Central Weather Bureau (
CWB
)of
Taiwan (Shin and Teng, 2001). It is currently the largest
well-recorded earthquake with
M
w
7.6. Four hundred and
forty-one strong-motion stations recorded the main event,
and 69 of those were at distances of less than 50 km from
the epicenter. We use three-component (north
south, east
west, and up
down) strong-motion records from the datset
collected by
CWB
. They classified the recorded accelero-
grams into four quality groups based on the existence of ab-
solute timing, preevents, and defects. For this analysis, we
use QA-class data (best for any studies), QB-class data (next
best, but no absolute timing), and a part of QC-class data
(covering the principal strong motions but not having pree-
vent or postevent data), which includes the preevent. The dis-
tribution of stations is shown in Figure 1.
Table 1 describes the crustal model for
P
-wave and
S
-wave velocity in central Taiwan (Ma
et al.
, 1996).
P
-wave
and
S
-wave arrival time for the predicted envelope are com-
puted with this
1D
layered crustal model. Because the origi-
nal seismic records reported incorrect universal time, we use
the data modified by Lee
et al.
(2001). They compared
picked
P
-wave arrival times with computed theoretical
P
-
wave arrival times. If the
P
-time residual was larger than
1 sec for accelerograms at the distance within 50 km, they
corrected the
P
-wave arrival time (Lee
et al.
, 2001). There-
fore, the error of the time stamp of the modified data is less
than 1 sec.
We use envelopes of acceleration records in horizontal
and vertical components. The horizontal components are cal-
culated by the square root of sum of squares of the east
west
and north
south components. Acceleration envelopes are ob-
tained by taking the maximum absolute amplitude of the
ground-motion time history over a 1-sec window. This
one second sampling acceleration envelopes reduces the
computational efforts.
Method
Ground-Motion Models
Cua and Heaton examined over 30,000 seismograms in
Southern California and developed relationships that predict
waveform envelopes (Cua, 2005). Each seismogram is para-
meterized with a set of 11 parameters, and each parameter is
characterized as a function of magnitude and epicentral dis-
tance. Therefore, the envelope ground-motion
E

t
j
M; R

at
time
t
can be expressed as a function of magnitude (
M
) and
epicentral distance (
R
). The detail of the ground-motion
models are explained in the Appendix. An observed ground-
motion envelope and the best
P
-wave,
S
-wave, and ambient
noise decomposition are shown in Figure 2. This ground-
motion model is similar to the traditional strong-motion at-
tenuation relationships (e.g., Boore
et al.
, 1993) with time
domain transition. However, unlike attenuation relationships
that have been developed primarily from large earthquakes,
the dataset of Cua and Heaton included earthquakes with
magnitudes ranging between 2 and 7. In order to characterize
this data, it was essential to utilize relationships in which the
distance decay of amplitudes change as a function of mag-
Table
1
P
- and
S
-Wave Velocity Model in Central Taiwan
*
Thickness (km)
V
P
(km
=
sec)
V
S
(km
=
sec)
1.0
3.50
2.00
3.0
3.78
2.20
5.0
5.04
3.03
4.0
5.71
3.26
4.0
6.05
3.47
8.0
6.44
3.72
5.0
6.83
3.99
Half-space
7.28
4.21
*
Ma
et al.
, 1996.
608
M. Yamada and T. Heaton
120 ̊00'
120 ̊15'
120 ̊30'
120 ̊45'
121 ̊00'
121 ̊15'
121 ̊30'
121 ̊45'
122 ̊00'
23 ̊15'
23 ̊30'
23 ̊45'
24 ̊00'
24 ̊15'
24 ̊30'
24 ̊45'
Shihtan fault
Chelungpu faul
t
0
50
km
Station
Epicenter
Subsource
0
1
2
3
4
5
6
7
−1
−2
−3
−4
C024
Taichung
Chiayi
Hualian
Figure
1.
The fault geometry and the station distribution of the Chi-Chi earthquake. The shaded area around the epicenter displays the
map projection of the fault geometry proposed by Ji
et al.
(2003). Small circles indicate the location of subsources determined from the
simulation in this article. The area within 50 and 100 km from the epicenter is shown by large circles. Stations used in this analysis are shown
by solid triangles. The polygon surrounding each station is the Voronoi cell for the station.
Figure
2.
Observed envelope for accelerogram and
P
-wave and
S
-wave envelopes for the ground-motion model defined in equation (A2)
(Cua, 2005).
Real-Time Estimation of Fault Rupture Extent Using Envelopes of Acceleration
609
nitude, similar to that which was introduced by Campbell
(1981). In particular, accelerations recorded close to a rup-
ture saturate at magnitudes larger than 6, whereas distant
sites do not demonstrate comparable saturation as a function
of magnitude.
Subsource Model
We introduce a
1D
source model to express the fault
finiteness based on the ground-motion models by Cua and
Heaton (Cua, 2005). The fault surface is divided into sub-
faults, and each subfault is represented by a single point
source, called subsources (Fig. 3). To simplify the problem,
we assume that the dimensions of all subsources are uniform.
Each source nucleates, and the
P
and
S
waves are radiated
when the rupture front arrives at the subsource.
The ground motion at a site is modeled as the combina-
tion of the responses of each subsource. An envelope of the
i
th subsource is expressed as follows:
E
i

t
j
M; R

E

t

t
i
j
M
i
;R
i

θ

;
(1)
where
M
i
is the magnitude of the subsource
i
,
R
i

θ

is the
distance between the station and subsource
i
, which is a func-
tion of the azimuthal angle (
θ
) of the rupture direction, and
t
i
is the time delay due to the rupture propagation.
For high-frequency motions with approximately random
phase, the square root of the sum of the squares of the en-
velope amplitudes from each subsource provides a good es-
timation of an acceleration envelope (Cocco and Boatwright,
1993):
E
total

t
j
M; R


X
N
1
i

0
E
2
i

t
j
M; R

X
N
2
i

1
E
2
i

t
j
M; R

v
u
u
t
;
(2)
where
N
1
and
N
2
are number of the point sources to the
north and the south, and
E
total

t
j
M; R

is the estimated en-
velope as a function of time.
E
i

t

is actually a fairly complex
function of time, magnitude, and distance, although its for-
ward calculation is very fast because it only involves analytic
functions (Cua, 2005; Cua and Heaton, 2007).
This model only works for high-frequency ground mo-
tions. Unlike longer-period ground motions, high-frequency
motions seem to be insensitive to either radiation pattern (Liu
and Helmberger, 1985) or directivity (Boatwright and Boore,
1982). Furthermore, near-source high-frequency motions
saturate as a function of magnitude. That is, near-source
high-frequency ground motions are independent of the am-
plitude of the slip for large earthquakes (Kanamori and Jen-
nings, 1978; Cua and Heaton, 2007).
Heaton and Hartzell (1989) pointed out that the assump-
tions of a Brune (1970) source spectrum combined with con-
stant stress drop leads to high-frequency energy radiated
from a subfault that is independent of the slip on the subfault.
A consequence of the fact that high-frequency near-source
ground motions can be modeled as random noise whose am-
plitude is independent of slip is that the high-frequency ra-
diated energy in earthquakes is proportional to the rupture
surface area. This is consistent with the observation of Boat-
wright (1982), who showed that high-frequency spectral ac-
celeration amplitudes are proportional to the root mean
square dynamic stress drop and the square root of the rupture
area. Our simple model for simulating high-frequency mo-
tions is also compatible with the observation of Hanks
and McGuire (1981): high-frequency ground accelerations
are remarkably similar from one event to another. Subsources
for our source model are evenly spaced, so the surface area
and high-frequency radiated energy corresponding to each
subsource are also constant. Based on this theoretical inter-
pretation, we estimated the ground-motion envelopes with
the subsource model for the 1999 Chi-Chi earthquake.
Figure 4 (top) shows an example of predicted envelopes
for vertical accelerations using the subsource model. It shows
the envelopes of the vertical acceleration record for each sub-
source with
M
6.0. Figure 4 (bottom) shows the time history
envelope of the accelerogram (vertical component) at the sta-
tion C024, a station on the foot wall side and 10 km from the
Chelungpu fault line (shown in Fig. 1, southwest of the epi-
center). Figure 4 also shows that the vertical acceleration en-
velopes predicted by the subsource model for the
VS-FS
(finite source) method fit the observed envelopes much better
station
station
fault
azimuthal
angle
N2 = 3
N1 = 4
epicenter
subsource
station
Figure
3.
Schematic diagram of the multiple source model. The
fault rupture is assumed to propagate from the epicenter at the con-
stant velocity
v
R
. The fault is parameterized by
θ
,
N
1
, and
N
2
,
where
θ
is the azimuthal angle of the fault and
N
1
and
N
2
are
the number of subsources north and south of the epicenter, respec-
tively. The ground motion at a station is expressed as a combination
of the envelope from each subsource.
610
M. Yamada and T. Heaton
than the envelopes predicted by the single source model for
the
VS-PS
(point source) method.
Even though the Ch-Chi rupture has large spatial varia-
tions in the amplitude of the slip, it appears that the high-
frequency accelerations can be modeled as a sum of the
radiation from a uniform tiling of the
M
6.0 subfaults, based
on the random-phase assumption and saturation with regard
to magnitude.
Comparison between Subsource Model
and Observations
Figures 5 and 6 are a comparison of observed envelopes
and predicted envelopes from the best-fit source model (ex-
plained in the next section). The model consists of 12 sub-
sources distributed along a line trending 17° clockwise from
north; there are seven subsources north of the epicenter and
four subsources to the south. The predicted acceleration en-
velopes for this model agree well with the observed enve-
lopes. Predicted envelopes of near-source stations have some
discrepancy depending on the source process, but predicted
envelopes of far-source stations fit the observation well.
The vertical predicted envelopes of the stations in the
epicentral region (e.g., stations T078, T079, T084, and
T089) are of particular interest. The predicted envelopes
overestimate these observed envelopes for the first 10 sec,
but then underestimates the observed records 20 sec after
the earthquake
s origin. The fact that the largest accelerations
in the epicentral region occurred 20 sec after the origin time
seems to indicate that there may have been some rupture
complexity in the hypocentral region; perhaps there was
an early aftershock in the epicentral region 20 sec after
the first rupture. Although this feature is noteworthy, it does
not have a significant effect on the inversions, because the
epicentral stations are less important for estimating azi-
muthal angle and length of the fault.
Note that there is a discrepancy between the predicted
and observed horizontal envelopes of the stations located
about 40 km north of the recognized northern terminus of
the Chelungpu fault rupture (e.g., stations T045, T047,
and T095). The observed records have a sharp pulse that ap-
pears about 40 sec after the event onset in the records of sta-
tions. Shin and Teng (2001) suggested that these large
accelerations were generated by a secondary rupture, perhaps
on the Shihtan fault.
Estimating the Fault Parameters
We performed an inversion of acceleration envelopes to
retrieve fault model parameters. We assume that the location
of the epicenter is already estimated from the
VS-PS
method,
and that the fault ruptures bilaterally from the epicenter with
constant rupture velocity. Thus, the time delay for each sub-
source rupture is the distance from the epicenter divided by
020406080100
0
50
100
150
acceleration (cm/s
2
)
Predicted (VS−FS)
020406080100
0
50
100
150
200
250
time (sec)
acceleration (cm/s
2
)
Observed
Predicted (VS−PS)
Predicted (VS−FS)
Figure
4.
Envelopes of vertical acceleration recorded at the station C024 for the Chi-Chi earthquake. Top: Predicted envelopes of the
vertical acceleration record for each subsource with
M
6.0. Bottom: Observed envelope (in dotted black line) and predicted envelopes of the
point source model in
VS-PS
method (in solid gray line) and of the multiple source model in
VS-FS
method (in solid black line).
Real-Time Estimation of Fault Rupture Extent Using Envelopes of Acceleration
611
the rupture velocity. Therefore, parameters that we need to
estimate from the observed data are the azimuthal angle
(
θ
) of the rupture direction,
N
1
and
N
2
, that are used to si-
mulate each of the segments of the bilateral rupture. Three
parameters are computed at each second using only the data
available at that time.
The best estimate of the model parameters minimizes the
residual sum of the squares (
RSS
) between observed ground-
motion envelopes and predicted envelopes from the sub-
source model. The misfit function as a measure of goodness
of fit is defined as follows:
RSS

t

X
ns
i

1
X
2
j

1
X
t
k

1

A
ijk

^
A
ijk

2
;
(3)
where
ns
is the number of stations,
t
is the time in 1-sec in-
tervals (
Δ
t

1
) from the event onset, and
A
ijk
and
^
A
ijk
are
observed and predicted envelopes of component
j
(horizon-
tal and vertical) at station
i
at time
k
Δ
t
.
This form of the misfit function tends to emphasize the
importance of fitting stations with large accelerations. That
is, distant stations have small observed and predicted accel-
erations, and even if there are serious discrepancies in the
ratio of the predicted and observed amplitudes, they will
have little impact on the inversion. In the course of this study,
we also tried inversions in which we defined misfit function
in terms of log of amplitudes. This misfit function empha-
sizes the ratio of predicted and observed amplitudes; large
amplitude data is no more important than small amplitude
data. However, we found that such a misfit function empha-
sized misfits in the coda for near-source data and furthermore
the distant data was often not well explained by our simple
descriptions of wave envelopes that have been developed to
explain the average effects of waves propagating through the
120 ̊00'
120 ̊15'
120 ̊30'
120 ̊45'
121 ̊00'
121 ̊15'
121 ̊30'
121 ̊45'
122 ̊00'
23 ̊00'
23 ̊15'
23 ̊30'
23 ̊45'
24 ̊00'
24 ̊15'
24 ̊30'
24 ̊45'
25 ̊00'
0
20
40
60
km
Taichung
Tainan
Chiayi
Hualian
Ilan
Station
Epicenter
Subsource
Predicted
Observed
60 secs.
C008
C012
C024
C034
C041
C052
C061
C074
C080
C082
C093
C099
C104
C107
C116
H019
H026
H033
H038
H041
H045
H056
H059
I063
I064
I067
K001
K054
N001
N014
N041
T026
T035
T038
T045
T059
T067
T068
T072
T076
T081
T084
T085
T088
T095
T111
T115
T119
T128
T147
Time (s) from event
Figure
5.
Predicted and observed envelopes in horizontal component. The red and black lines are the predicted and observed envelopes,
respectively. The locations of the subsources estimated from the best fit model are shown in a small yellow circles. The area within 50 and
100 km from the epicenter are shown by large circles. Only characteristic records of the stations are shown in this figure.
612
M. Yamada and T. Heaton
crust. That is, it is important to emphasize the data from the
near-source records and a logarithmic misfit function was not
appropriate to recover the timing and location of the rupture.
Our parameterization scheme has the advantage that
we characterize the source with relatively few parameters

θ
;N
1
;N
2

, none of which require high precision estimates.
However, for this strategy to be effective, we will need to
solve a nonlinear inverse problem in real time. In this study,
we solve this nonlinear inverse problem by using the neigh-
borhood algorithm (Sambridge, 1999a,b). We recognize that
other inverse techniques may ultimately be chosen for real-
time applications. However, because the purpose of this
study is to determine the effectiveness of our parameteriza-
tion, we use the neighborhood algorithm to characterize and
solve this nonlinear inverse problem.
The neighborhood algorithm is a direct search method
for finding models of acceptable data fit in a multidimen-
sional parameter space. We generate samples in the para-
meter space and draw the Voronoi cells for these samples.
Voronoi cells are nearest-neighbor regions defined under a
suitable distance norm, and the shape and the size of each
Voronoi cell is determined by the sample distribution in
the parameter space. See Figure 1 as an example of Voronoi
cell as used to define the nearest neighbors to seismic sta-
tions. We calculate the misfit function for each sample
and choose the model with the minimum misfit. New sam-
ples are generated by performing a uniform random walk in
the chosen Voronoi cell. By repeating these steps, we will
find a set of samples that identifies those regions of the para-
meter space that provide the best fit to the data. This is an
120 ̊00'
120 ̊15'
120 ̊30'
120 ̊45'
121 ̊00'
121 ̊15'
121 ̊30'
121 ̊45'
122 ̊00'
23 ̊00'
23 ̊15'
23 ̊30'
23 ̊45'
24 ̊00'
24 ̊15'
24 ̊30'
24 ̊45'
25 ̊00'
0
20
40
60
km
Taichung
Tainan
Chiayi
Hualian
Ilan
Station
Epicenter
Subsource
Predicted
Observed
60 secs.
C012
C024
C034
C041
C052
C061
C074
C080
C082
C093
C099
C104
C107
C116
H019
H026
H033
H038
H041
H045
H056
H059
I063
I064
I067
K001
K054
N001
N014
N041
T026
T035
T038
T045
T059
T067
T068
T072
T076
T081
T084
T085
T088
T095
T111
T115
T119
T128
T147
Time (s) from event
Figure
6.
Predicted and observed envelopes in vertical component. The red and black lines are the predicted and observed envelopes,
respectively. The locations of the subsources estimated from the best -fit model are shown in a small yellow circles. The areas within 50 and
100 km from the epicenter are shown by large circles. Only characteristic records of the stations are shown in this figure.
Real-Time Estimation of Fault Rupture Extent Using Envelopes of Acceleration
613
approach for constructing the posterior probability density
function from the ensemble samples based on the Voronoi
cell concept (Yamada, 2007).
Results of the Simulation for the Chi-Chi Earthquake
We have run many different inversions by varying both
the inversion parameters and the data sampling, such as the
number of records used for the inversion, the components of
the records, etc. In this section, we show the results of the
main simulation. We use the horizontal and vertical records
of the stations within 120 km of the epicenter. To simplify the
problem, we assume each subsource has the same magnitude
of 6.0 and the same spacing of 10 km.
Result of the Fault Parameter Estimate
Figure 7 shows the estimation results for three para-
meters: azimuthal angle of fault line (
θ
), number of the point
sources to the north (
N
1
), and number of point sources to the
south (
N
2
). These three parameters are computed at each
second using only the data available at that time. The esti-
mation is updated every second as the ground-motion data
are observed. Although it does a good job at characterizing
the rupture length and timing, we see that it is difficult to
resolve
θ
until 15 sec after the event onset because the event
can be approximated as a point source at the beginning. The
estimated
θ
at 15 sec is about 15°, and it increases gradually
after 20 sec due to an impulsive acceleration arrival at station
C080, which is located at the south of the epicenter. Esti-
mates of
θ
stabilize at about 14° with respect to additional
data after 27 sec. There is an additional small shift at
45 sec, at which point the inversion achieves its final solution
of 17°, which compares favorably with the observed average
fault strike of the Chelungpu fault rupture (the fault trace is
shown in Fig. 1).
Because the subsources are equally spaced, the length of
the fault is represented by the number of the point sources to
the north (
N
1
) and to the south (
N
2
). Figure 7 (bottom)
shows values of
N
1
and
N
2
as a function of time after
the origin. From the figure, we can see the fault length grows
bilaterally along the dashed black lines. Around 30 sec, the
rupture stops growing to the south. It also stops to the north
temporarily, but it grows again around 40 sec. This is due to
the delayed high-frequency radiation at stations north of the
Chenlungpu surface rupture and may have been caused by
rupture on the Shihtan fault. Even though the result of the
simulation fits the actual location of the fault accurately,
the subsource model does not consider rupture jumping dis-
locations (i.e., the rupture at the adjacent active faults trig-
gered by the mainshock) (Shin and Teng, 2001). The final
result shows seven point sources to the north and four point
sources to the south. This fault length is comparable to the
total length from the Chelungpu fault to the Shihtan fault in
Figure 1.
The best-fit source model consists of 12 subsources dis-
tributed along a line trending 17° clockwise from north; there
are seven subsources north of the epicenter and four sub-
sources to the south. That is, the best-fitting source parameter
is given by (
θ

17
°,
N
1

7
, and
N
2

4
). The predicted
acceleration envelopes for this model agree well with the ob-
served envelopes, as shown in Figures 5 and 6. Predicted en-
velopes of near-source stations have some discrepancies
depending on the source process, but predicted envelopes
of far-source stations fit the observation well.
Sensitivity Analysis of Station Sampling
and Rupture Velocity
Sensitivity of the station distribution and rupture velo-
city are also tested. We assume a constant rupture velocity
of
2
:
0
km
=
sec to construct the predicted envelopes from sub-
sources. In order to check the sensitivity of the parameter
estimate to the rupture velocity, we run four simulations with
different rupture velocities. Figure 8 shows the estimated
parameters,
N
1
and
N
2
, for rupture velocities varying from
2.0 to
3
:
5
km
=
sec. The simulation results are relatively in-
sensitive to the choice of rupture velocity. Because the most
probable rupture velocity minimizes the misfit between ob-
served envelopes and predicted envelopes, we can include
this rupture velocity as a parameter to be estimated for future
analysis.
In order to test the effect of station distribution, some of
the stations are removed from the dataset. By randomly sam-
pling from the total of 239 stations, we created two subsets of
the dataset: subset A consists of 126 records with an even
station code number, and subset B includes 56 records with
a station code ending in 6 or 8 (e.g., T078). Although the
station distribution is not homogeneous, the average station
density is
214
km
2
=
station for subset A, and
482
km
2
=
station for subset B. The stations are located in an area of
about
27
;
000
km
2
. Even though the station density is differ-
ent, the estimated parameters are quite similar. In Figure 7,
the time series of
θ
and
N
2
for full dataset, subset A, and
subset B are almost the same.
N
1
for subsets A and B stays
around 5 after 30 sec, and the increase observed in full data-
set due to the Shihtan fault rupture does not appear. The rea-
son is that several near-source stations of the Shihtan fault
have an odd-numbered station code and are not included
in this analysis (e.g., T045, T047, and T095). Considering
that the rupture of the Shihtan fault is quite small compared
to that of the Chelungpu fault, Subsets A and B can express
the Chi-Chi earthquake rupture well. The
VS-FS
method for
large earthquakes works well even if the station density is
reduced to a quarter of the original density as long as the
station distribution is uniform.
Geometry of the Parameter Space
We have solved the optimization problem in parameter
space (
θ
,
N
1
, and
N
2
) by the neighborhood algorithm. Here,
we discuss the geometry of the parameter space.
614
M. Yamada and T. Heaton
Figure 9 shows the error surface of
θ
and
N
1
at a fixed
N
2
of 4 and assuming that all data is used in the inversion.
The surface is smooth and has a deep and narrow valley at
θ

10
°. The solution easily converges to this minimum.
Figure 10 shows the error surface of
N
1
and
N
2
at a fixed
θ
of 17°. The surface is very smooth in both
N
1
and
N
2
directions. The global minimum is very sensitive to the
choice of the dataset, as shown in the results of subsets A
and B.
Contour maps of the error surface of
N
1
and
N
2
at 10-
sec intervals are shown in Figure 11.
θ
is fixed at 17°, which
is the optimal final solution. At 10 sec, the minimum of this
error surface is

N
1
;N
2

0
;
1

. However, it is not the glo-
bal minimum in the parameter space because
θ

17
° is not
the optimal solution at 10 sec. At 20 and 30 sec, the mini-
mum of the error surface is at the maximum
N
1
and
N
2
in the
possible parameter space, even though
θ
is not optimal.
There is a strong possibility that the rupture is still ongoing
at this point. At 40 sec, the minimum of the contour is around

N
1
;N
2

5
;
4

and it suggests that the rupture has
stopped rupturing toward the south. After 40 sec, the shape
of contour map does not change much, and the elliptic shape
0
10
20
30
40
50
60
50
0
50
100
time(s)
azimuthal angle
θ
(deg)
full dataset (239 stations)
subset A (126 stations)
subset B (56 stations)
0
10
20
30
40
50
60
15
10
5
0
5
10
15
time(s)
(south) N2 N1 (north)
V
r
=2.0 km/s
Number of the subsources
full dataset (239 stations)
subset A (126 stations)
subset B (56 stations)
Figure
7.
Time series of the estimated parameters,
θ
,
N
1
, and
N
2
for each dataset. Time is relative to the origin. The parameters are
computed at each second using only the data available at that time. The broken lines are the best estimates based on the fault model proposed
by Ji
et al.
(2003). Top: Time series estimations for
θ
. Bottom: Time series estimations for
N
1
and
N
2
. The solid thin lines are the upper limits
for
N
1
and
N
2
for the rupture velocity
2
km
=
sec.
Real-Time Estimation of Fault Rupture Extent Using Envelopes of Acceleration
615
of the smallest contour indicates that
N
2
is determined un-
iquely, but that considerable uncertainty about
N
1
remains.
Discussion and Conclusions
We outlined a strategy to estimate slip in time and space
for an ongoing earthquake rupture. A key aspect of this strat-
egy is to map the location of the rupture using envelopes of
high-frequency acceleration data. Once the location of the
rupture is estimated, long-period displacement data can be
projected back onto the fault to determine the slip in real time
(see Yamada, 2007).
Our strategy for using high-frequency radiation to deter-
mine the timing and length of the rupture relies on the ob-
servation that high-frequency seismic waves can be modeled
as random-phase waves whose total radiated energy scales
linearly with the rupture area. By using this assumption,
we show that we can simulate the ground motion of a large
earthquake by tiling the surface of the large event with smal-
ler events and then summing the random-phase signals from
the smaller events. In our example of the Chi-Chi earthquake,
0
10
20
30
40
50
60
100
50
0
50
100
time(s)
azimuthal angle
θ
(deg)
Vr = 2.0km/s
Vr = 2.5km/s
Vr = 3.0km/s
Vr = 3.5km/s
0
10
20
30
40
50
60
15
10
5
0
5
10
15
time(s)
(south) N2 N1 (north)
V
r
=2.0 km/s
V
r
=3.5 km/s
Number of the subsources
Vr = 2.0km/s
Vr = 2.5km/s
Vr = 3.0km/s
Vr = 3.5km/s
Figure
8.
The estimated parameters,
N
1
and
N
2
for different rupture velocities. The solid thin lines are the upper limits for
N
1
and
N
2
for
the rupture velocity 2 and
3
:
5
km
=
sec. The broken lines are the best estimates based on the fault model proposed by Ji
et al.
(2003). Time is
relative to the origin. The parameters are computed at each second using only the data available at that time.
616
M. Yamada and T. Heaton
0
2
4
6
8
10
90
60
30
0
30
60
90
40
45
50
55
60
θ
(deg)
N1
misfit
45
50
55
Figure
9.
Error surface of
θ
and
N
1
at the fixed
N
2

4
at
60 sec after the origin time. Because the surface is peaked around
θ

0
, it is easy to converge in
θ
. However, the optimal
N
1
will
change easily depending on the misfit function (see equation 3).
0
2
4
6
8
10
0
2
4
6
8
10
40
45
50
55
60
65
N1
N2
misfit
45
50
55
60
Figure
10.
Error surface of
N
1
and
N
2
at the fixed
θ

17
at
60 sec after the origin time. Because the surface is smooth in both
N
1
and
N
2
directions; the optimal solution is sensitive to a small
disturbance.
N1
N2
T = 10sec
0
2
4
6
8
10
0
2
4
6
8
10
N1
N2
T = 20sec
0
2
4
6
8
10
0
2
4
6
8
10
N1
N2
T = 30sec
0
2
4
6
8
10
0
2
4
6
8
10
N1
N2
T = 40sec
0
2
4
6
8
10
0
2
4
6
8
10
N1
N2
T = 50sec
0
2
4
6
8
10
0
2
4
6
8
10
N1
N2
T = 60sec
0
2
4
6
8
10
0
2
4
6
8
10
Figure
11.
Contour maps of the error surface of
N
1
and
N
2
at the fixed
θ

17
. The maps are shown in 10-sec intervals. The blank area
in the boxes is the region where there is no solution due to the constraint that the rupture velocity is less than
2
km
=
sec.
Real-Time Estimation of Fault Rupture Extent Using Envelopes of Acceleration
617