M. Yamada and T. Heaton

1

Real-time Estimation of Fault Rupture Extent using
Envelopes of Acceleration
Masumi Yamada1, Thomas Heaton2
masumi@eqh.dpri.kyoto-u.ac.jp, heaton t@caltech.edu

Abstract:
We present a new strategy to estimate the geometry of a rupture on a finite fault in realtime 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 relationships. This paper extends the VS method
to large earthquakes where fault finiteness is important. We propose a new model to simulate
high-frequency motions from earthquakes with large rupture dimension: the envelope of highfrequency ground motion from a large earthquake can be expressed as a root-mean-squared
combination of envelope functions from smaller earthquakes. We use simulated envelopes ofground 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 parameters
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 tothe minimum uniquely, and the estimation agrees with the actual Chi-Chi earthquake fault
geometry quite well. The rupture direction can be estimated at 10 seconds after the event
onset and the final solution is achieved after 20 seconds. While this methodology seems quite
promising for warning systems, it only works well when there is an adequate distribution ofnear-source stations.

1

Kyoto University,
Gokasyo, Uji, 611-0011, Japan2
California Institute of Technology,
MC104-44, 1200 E. California Blvd. Pasadena CA 91125 USA

M. Yamada and T. Heaton

12Introduction

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 various
early warning methods have been proposed from seismologists and engineers (Nakamura,
1988; Allen and Kanamori, 2003; Odaka et al., 2003; Wu and Kanamori, 2005). Currently,
the most ambitious system is the earthquake early warning system provided by the Japan
Meteorological Agency, which is released to the public in October, 2007. The news of thesystem was broadcasted widely and 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, anearthquake early warning system must collect and quickly analyze seismic data in a mannerthat can be used to predict future shaking. In principle, this could be achieved by using
the present value of an approximately known wavefield as a boundary condition to predict
future wavefields using Navier’s equation (Baker et al., 2005). However, from a practical
viewpoint, there are advantages to data analysis schemes that involve characterization of theearthquake source. Predictions of future shaking can be achieved by utilizing the extensive
existing work on predicting ground shaking from seismic sources. Ideally, an early warning
system would provide the best estimate of slip in time and space that can be deduced fromseismic data available at any given instant in time. We introduce a two-step strategy toaccomplish this; 1) we determine the spatial and temporal extent of an ongoing ruptureby analyzing waveform envelopes of high-frequency shaking, 2) we determine approximate
slip from simple projections of long-period shaking onto the approximately known locationof the rupture. It seems likely that true ground displacement may be very helpful in thereal-time determination of fault slip. Since it is difficult to obtain displacement from inertial
seismometers (Clinton, 2004), the future incorporation of real-time high-sample-rate GPS
into early warning systems may be quite important. A more complete description of ourstrategy to determine distributions in real time is given by Yamada (2007).
Cua and Heaton developed the Virtual Seismologist (VS) method (Cua, 2005). It isa Bayesian approach to seismic early warning designed for modern seismic networks, andproposed for small to moderate earthquakes with ruptures that can be approximately modeled as a point source. The VS algorithm uses an envelope attenuation relationship andthe 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 magnitude and location (Cua, 2005).
It gives the best estimate of an earthquake property in terms of a probability density function. The Bayesian approach is a scheme to emulate human capabilities to judge complex
information by modeling uncertainty 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 realtime (Nakamura, 2001), most of earthquake early warning systems focus on estimating epicenters and magnitudes of earthquakes, not the fault geometry (Allen and Kanamori, 2003).
However, for large earthquakes, rupture length can be on the order of tens to hundreds ofkilometers, and the inhomogeneous slip distribution significantly affects the ground motion

M. Yamada and T. Heaton

3

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 fromthe epicenter when using a scheme that only characterizes the earthquake as a point source.
In this paper, we introduce a methodology to determine the direction and length of
the rupture from high-frequency ground motions. We construct a model to describe highfrequency motions from earthquakes with large fault dimension. By applying this model tohigh-frequency ground motions in real time, it becomes possible to estimate the geometry ofan ongoing earthquake.

2

Data

The data for this analysis is the strong motion dataset from the 21 September 1999 Chi-Chi
Earthquake that occurred in central Taiwan (Lee et al., 2001). The epicenter was locatedat 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 earthquakewith moment magnitude 7.6. 441 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 (NS,EW, and UD) strong motion records from the data set collected by CWB. They classified therecorded accelerograms into four quality groups based on the existence of absolute timing,
pre-events, and defects. For this analysis, we use QA-class data (best for any studies), QBclass data (next best but no absolute timing) and a part of QC-class data (covering theprincipal strong motions but not having pre-event or post-event data) which includes thepre-event. The distribution of stations is shown in Figure 1.
Table 1 describes the crustal model for P-wave and S-wave velocity in central Taiwan (Maet al., 1996). P-wave and S-wave arrival time for the predicted envelope are computed withthis 1-D layered crustal model. Since the original 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 second for accelerograms at the distance within 50 km, they corrected the P-wave arrival
time (Lee et al., 2001). Therefore, the error of the time stamp of the modified data is lessthan 1 second.
We use envelopes of acceleration records in horizontal and vertical components. The
horizontal components are calculated by the square root of sum of squares of the EW and
NS components. Acceleration envelopes are obtained by taking the maximum absolute
amplitude of the ground motion time history over a 1 second window. This one second
sampling acceleration envelopes reduces the computational efforts.

4

M. Yamada and T. Heaton

3

Method3.1

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 parameterized with a set of eleven parameters, and each parameter is characterised as a function ofmagnitude and epicentral distance. Therefore, the envelope ground motion E(t|M, R) attime 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 groundmotion envelope and the best P-wave, S-wave, and ambient noise decomposition are shownin Figure 2. This ground motion model is similar to the traditional strong motion attenuation relationships (e.g., Boore et al., 1993) with time domain transition. However, unlike
attenuation relationships that have been developed primarily from large earthquakes, thedata set of Cua and Heaton included earthquakes with magnitudes ranging between 2 and7. In order to characterize this data, it was essential to utilize relationships in which the
distance decay of amplitudes change as a function of magnitude, similar to that which wasintroduced by Campbell (1981). In particular, accelerations recorded close to a rupture
saturate at magnitudes larger than 6, whereas distant sites do not demonstrate comparable
saturation as a function of magnitude.

3.2

Subsource Model

We introduce a one-dimensional 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
subfaults, and each subfault is represented by a single point source, called ”subsources”(Figure 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 combination of the responses of eachsubsource. An envelope of the ith subsource is expressed as follows:
Ei (t|M, R) = E(t − ti |Mi , Ri (θ))(1)

where Mi is the magnitude of the subsource i, Ri (θ) is the distance between the station andsubsource i which is a function of the azimuthal angle (θ) of the rupture direction, and ti isthe time delay due to the rupture propagation.
For high-frequency motions with approximately random phase, the square root of the sumof the squares of the envelope amplitudes from each subsource provides a good estimationof an acceleration envelope (Cocco and Boatwright, 1993):
v
u N1
N2
uX
X
t
2
Etotal (t|M, R) =Ei2 (t|M, R)(2)
Ei (t|M, R) +
i=0

i=1

M. Yamada and T. Heaton

5

where N1 and N2 are number of the point sources to the north and the south and Etotal (t|M, R)
is the estimated envelope as a function of time. Ei (t) is actually a fairly complex functionof time, magnitude, and distance, although its forward calculation is very fast since it only
involves analytic functions (Cua, 2005; Cua and Heaton, 2006).
This model only works for high-frequency ground motions. Unlike longer-period groundmotions, high-frequency motions seem to be insensitive to either radiation pattern (Liu andHelmberger, 1985) or directivity (Boatwright and Boore, 1982). Furthermore, near-sourcehigh-frequency motions saturate as a function of magnitude. That is, near-source highfrequency ground motions are independent of the amplitude of the slip for large earthquakes(Kanamori and Jennings, 1978; Cua and Heaton, 2006).
Heaton and Hartzell (1989) pointed out that the assumptions of a Brune (1970) source
spectrum combined with constant 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 highfrequency near-source ground motions can be modeled as random noise whose amplitude isindependent of slip is that the high-frequency radiated energy in earthquakes is proportional
to the rupture surface area. This is consistent with the observation of Boatwright (1982)
who showed that high-frequency spectral acceleration amplitudes are proportional to theroot-mean-square (rms) dynamic stress drop and the square root of the rupture area. Our
simple model for simulating high-frequency motions is also compatible with the observation ofHanks and McGuire (1981) that high-frequency ground accelerations are remarkably similarfrom 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 interpretation, we estimated the ground motion envelopes with thesubsource 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 eachsubsource with magnitude 6.0. Figure 4 (bottom) shows the time history envelope of theaccelerogram (vertical component) at the station C024, a station on the foot wall side and
10 km from the Chelungpu fault line (shown in Figure 1, southwest of the epicenter). Figure
4 also shows that the vertical acceleration envelopes predicted by the subsource model for
the VS-FS (finite source) method fit the observed envelopes much better 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 variations in the amplitude of the slip,
it appears that the high-frequency accelerations can be modeled as a sum of the radiationfrom a uniform tiling of the magnitude 6.0 subfaults, based on the random-phase assumption
and saturation with regard to magnitude.

3.3

Comparison between Subsource Model and Observations

Figures 5 and 6 are a comparison of observed envelopes and predicted envelopes from thebest fit source model (explained in the next section). The model consists of 12 subsources
distributed along a line trending 17 degrees clockwise from north; there are 7 subsources
north of the epicenter and 4 subsources to the south. The predicted acceleration envelopes

6

M. Yamada and T. Heaton

for this model agree well with the observed envelopes. 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 seconds, but then underestimates the observed
records 20 seconds after the earthquakes origin. The fact that the largest accelerations in
the epicentral region occurred 20 seconds 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 seconds after the first rupture. Although thisfeature is noteworthy, it does not have a significant effect on the inversions since the epicentral
stations are less important for estimating azimuthal 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 appears about 40 seconds after the event onset in the records ofstations. Shin and Teng (2001) suggested that these large accelerations were generated by asecondary rupture, perhaps on the Shihtan fault.

3.4

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, andthat the fault ruptures bilaterally from the epicenter with constant rupture velocity. Thus,
the time delay for each subsource rupture is the distance from the epicenter divided by therupture velocity. Therefore, parameters that we need to estimate from the observed data are
the azimuthal angle (θ) of the rupture direction, N1 and N2, that are used to simulate eachof the segments of the bilateral rupture. Three parameters are computed at each secondusing 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 subsource model. The misfit function as a measure of goodness of fit is defined as follows:
RSS(t) =

ns X
2 X
t
X

(Aijk − Âijk )2(3)

i=1 j=1 k=1

where ns is the number of stations, t is the time in 1 second intervals (∆t = 1) from the event
onset, and Aijk and Âijk are observed and predicted envelopes of component j (horizontal
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 accelerations 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

M. Yamada and T. Heaton

7function emphasizes 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 emphasized 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 beendeveloped to explain the ”average” effects of waves propagating through the 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 (θ, N1, N2), 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 Neighborhood Algorithm(Sambridge, 1999a; Sambridge, 1999b). We recognize that other inverse techniques may
ultimately be chosen for real-time applications. However, since the purpose of this study is
to determine the effectiveness of our parameterization, 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 acceptabledata fit in a multidimensional parameter space. We generate samples in the parameter spaceand 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 isdetermined 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 stations. We calculate
the misfit function for each sample and choose the model with the minimum misfit. New
samples 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 theparameter space that provide the best fit to the data. This is an approach for constructing
the posterior probability density function from the ensemble samples based on the Voronoi
cell concept (Yamada, 2007).

4

Results of the Simulation for the Chi-Chi Earthquake

We have run many different inversions by varying both the inversion parameters and thedata 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 thehorizontal 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 spacingof 10km.

4.1

Result of the Fault Parameter Estimate

Figure 7 shows the estimation results for three parameters (azimuthal angle of fault line(θ), number of the point sources to the north (N1) and to the south (N2)). These threeM. Yamada and T. Heaton

8

parameters are computed at each second using only the data available at that time. The
estimation 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 seconds after the event onset since the event can be approximated asa point source at the beginning. The estimated θ at 15 seconds is about -15 degrees and it
increases gradually after 20 seconds due to a impulsive acceleration arrival at station C080
which is located at the south of the epicenter. Estimates of θ stabilize at about 14 degrees
with respect to additional data after 27 seconds. There is an additional small shift at 45
seconds, at which point the inversion achieves its final solution of 17 degrees, which compares
favorably with the observed average fault strike of the Chelungpu fault rupture (the fault
trace is shown in Figure 1).
Since the subsources are equally spaced, the length of the fault is represented by the
number of the point sources to the north (N1) and to the south (N2). Figure 7 (bottom)
shows values of N1 and N2 as a function of time after the origin. From the figure, we cansee the fault length grows bilaterally along the dashed black lines. Around 30 seconds, therupture stops growing to the south. It also stops to the north temporarily, but it grows again
around 40 seconds. 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, thesubsource model does not consider ”rupture jumping dislocations” (i.e. the rupture at theadjacent active faults triggered by the main shock) (Shin and Teng, 2001). The final resultshows 7 point sources to the north and 4 point sources to the south. This fault length iscomparable to the total length from the Chelungpu fault to the Shihtan fault in Figure 1.
The best fit source model consists of 12 subsources distributed along a line trending 17
degrees clockwise from north; there are 7 subsources north of the epicenter and 4 subsources
to the south. That is, the best fitting source parameter is given by (θ=17 degrees, N1=7,
N2=4). The predicted acceleration envelopes for this model agree well with the observed
envelopes, as shown in Figures 5 and 6. Predicted envelopes of near-source stations havesome discrepancies depending on the source process, but predicted envelopes of far-source
stations fit the observation well.

4.2

Sensitivity Analysis of Station Sampling and Rupture Velocity

Sensitivity of the station distribution and rupture velocity are also tested. We assume aconstant rupture velocity of 2.0 km/s to construct the predicted envelopes from subsources.
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,
N1 and N2, for rupture velocities varying from 2.0 km/s to 3.5 km/s. The simulation results
are relatively insensitive to the choice of rupture velocity. Since the most probable rupture
velocity minimizes the misfit between observed 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 fromthe dataset. By randomly sampling from the total of 239 stations, we created two subsetsof the dataset: subset A consists of 126 records with an even station code number andM. Yamada and T. Heaton

9subset B includes 56 records with a station code ending in 6 or 8 (e.g. T078). Although thestation distribution is not homogeneous, the average station density is 214 km2 /station forthe subset A, and 482 km2 /station for the subset B. The stations are located in an area of
about 27000 km2 . Even though the station density is different, the estimated parameters
are quite similar. In Figure 7, the time series of θ and N2 for full dataset, subset A, andsubset B are almost the same. N1 for subsets A and B stays around 5 after 30 seconds, andthe increase observed in full dataset due to the Shihtan fault rupture does not appear. The
reason is that several near-source stations of the Shihtan fault have an odd number station
code and are not included in this analysis (e.g. T045, T047, and T095). Considering thatthe 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 forlarge earthquakes works well even if the station density is reduced to a quarter of the originaldensity, as long as the station distribution is uniform.

4.3

Geometry of the Parameter SpaceWe have solved the optimization problem in parameter space (θ, N1, and N2) by the Neighborhood Algorithm. Here, we discuss the geometry of the parameter space.
Figure 9 shows the error surface of θ and N1 at a fixed N2 of 4 and assuming that alldata is used in the inversion. The surface is smooth and has a deep and narrow valley atθ = 10 degrees. The solution easily converges to this minimum. Figure 10 shows the error
surface of N1 and N2 at a fixed θ of 17 degrees. The surface is very smooth in both N1 andN2 directions. The global minimum is very sensitive to the choice of the dataset, as shownin the results of subsets A and B.
Contour maps of the error surface of N1 and N2 at 10 second intervals are shown in
Figure 11. θ is fixed at 17 degrees which is the optimal final solution. At 10 seconds, theminimum of this error surface is (N1, N2) = (0, 1). However, it is not the global minimumin the parameter space since θ = 17 degrees is not the optimal solution at 10 seconds. At
20 and 30 seconds, the minimum of the error surface is at the maximum N1 and N2 in thepossible parameter space even though θ is not optimal. There is a strong possibility thatthe rupture is still ongoing at this point. At 40 seconds, the minimum of the contour isaround (N1, N2) = (5, 4) and it suggests that the rupture has stopped rupturing toward thesouth. After 40 seconds, the shape of contour map does not change much, and the elliptic
shape of the smallest contour indicates that N2 is determined uniquely, but that considerable
uncertainty about N1 remains.

5

Discussion and Conclusions

We outlined a strategy to estimate slip in time and space for an ongoing earthquake rupture.
A key aspect of this strategy is to map the location of the rupture using envelopes of highfrequency 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).

10

M. Yamada and T. Heaton

Our strategy for using high-frequency radiation to determine the timing and length of
the rupture relies on the observation 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 earthquakeby tiling the surface of the large event with smaller events and then summing the random
phase signals from the smaller events. In our example of the Chi-Chi earthquake, we showedthat a sum of 10 km interval magnitude 6.0 subevents provided a good prediction of theacceleration envelopes for this earthquake. In order to turn this simulation into a real-timeinverse, we parameterize the rupture with as a linear alignment of magnitude 6.0 earthquakes.
We then invert for the azimuth angle of the alignment as well as two integers N1 and N2,
which are the number of additional 10 km patches in the positive and negative directions
from the epicenter, respectively.
The best estimate of the model parameters minimizes the residual sum of the squares between observed ground motion envelopes and predicted envelopes from the subsource model(in equation (3)). This misfit function with linear amplitudes of ground motions can provide
better estimates than that of logarithmic amplitudes, since the linear misfit function tends
to emphasize the importance of fitting stations with large amplitudes.
Our study of the Chi-Chi data set indicates that it is more difficult to determine rupture
length than it is to determine rupture azimuth. Furthermore, for this method to work well,an adequate near-source station distribution is important.
We are currently investigating alternate strategies to determine the temporal evolutionof the rupture. In particular, we have developed a method to distinguish whether or not aparticular seismic record is from a station within 10 km of a rupture surface (Yamada et al.,,2007). Real-time mapping of an going rupture using this strategy becomes a simple matter oftracking the spatial evolution of near-source seismic stations. Although this strategy appears
promising, it requires adequate station coverage to track near-source stations.

Appendix: Ground Motion Models
Cua and Heaton examined over 30,000 seismograms in Southern California and developed
relationships that predict waveform envelopes as a function of magnitude, distance andstation corrections (Cua, 2005). First, the envelopes of the ground motions are modeled as
a combination of the envelopes of P-wave, S-wave, and ambient noise.
q2
Eobserved (t) = EP2 (t) + ES2 (t) + Eambient
+ǫ
(4)
where
Eobserved (t)EP (t)ES (t)
Eambient
ǫ

=
=
=
=
=

envelope of observed ground motion
envelope of P-wave
envelope of S-wave and later-arriving phases
ambient noise at the site
difference between predicted and observed envelope

11

M. Yamada and T. Heaton

The ambient noise, Eambient , for a given time history is modeled as a station constant. The
P- and S-wave envelopes, EP (t) and ES (t), are defined by a rise time (triseP and triseS ), a
constant amplitude (AP and AS ), a duration (∆tP and ∆tS ), and two decay parameters (γPand γS ) and (τP and τS ) respectively. The general form of the envelope function is:0; t < Ti



A
i,j

(t − Ti ); Ti ≤ t < Ti + trisei,jtrisei,j
Ei,j (t) =
(5)

Ai,j
; Ti + trisei,j ≤ t < Ti + trisei,j + ∆ti,j
1
 Ai,j; t ≥ Ti + trisei,j + ∆ti,j(t−Ti −trise −∆ti,j +τi,j )γi,j
i,j

where
i = P-, S-wave
Ti = P-, S-wave arrival times
j = acceleration, velocity, and displacement(horizontal and vertical components)

They parameterized each seismogram as a set of eleven parameters (five for the P-wave
envelope, five for the S-wave envelope, and one for the ambient noise). Furthermore, eachparameter is described by magnitude, distance, log of distance, and site dependent constants
based on the traditional strong motion attenuation relationships (Campbell, 1981; Boore etal., 1993). The functional forms which describe the P- and S- wave envelope functions are
given below:

log10 Aij =aij M + bij (R1 + Cij (M)) + dij log10 (R1 + Cij (M)) + eij + ǫij ,log10 Bij =aij M + bij R1 + dij log10 R1 + eij + ǫij ,

(6)
(7)

where
i =P-, S-wave,
j =horizontal and vertical acceleration, velocity, and displacement,
Aij =ground motion envelope amplitude,
Bij =rise time (trise ), duration (∆T ), and decay parameters (τ , γ),
M =local magnitude (Mw for M > 5.0),
R = epicentral distance in km for M < 5,closest distance to fault for M > 5.0 (when available),p
R1 = (R2 + 9),
Cij (M) =(arctan(M − 5) + 1.4)(c1ij exp(c2ij (M − 5))),aij , bij , c1ij , c2ij , dij , eij =regression constants,
ǫij =statistical (or prediction) error, ∼ NID(0, σ 2 ).

M. Yamada and T. Heaton

12

The Aij s are the ground motion envelope amplitudes (P- or S-wave) from fitting equations 4 and 5 to the observed ground motion envelopes in the database. The Bij s are theparameters characterizing the envelope function (trise , ∆T , τ , and γ). Coefficients in equations 6 and 7 are determined by regression analysis of the database using the Neighborhood
Algorithm (Cua, 2005).

Acknowledgments
The authors appreciate Dr. Georgia Cua for providing the codes and helpful advice for
the Virtual Seismologist method and the ground motion attenuation relationships. We also
thank Dr. Chen Ji for providing the fault model for the 1999 Chi-Chi earthquake. Part of the
figures are generated by Generic Mapping Tools (Wessel and Smith, 1991). We appreciate
the CWB of Taiwan and Dr. Yih-Min Wu for providing the ground-motion data. This
research was supported by the U.S. Geological Survey.

References- Allen, R. M., and H. Kanamori (2003). The potential for earthquake early warning inSouthern California, Science, 300, 685-648.
- Baker, T., R. Granat, and R. W. Clayton (2005). Real-time Earthquake Location
Using Kirchhoff Reconstruction, Bull. Seism. Soc. Am. 75, 699-707.
- Boatwright, J. (1982). A dynamic model for far-field acceleration, Bull. Seism. Soc.
Am. 72, 1049-1068.
- Boatwright, J., and D. M. Boore (1982). Analysis of the ground accelerations radiated by the 1980 Livermore Valley earthquakes for directivity and dynamic source
characteristics, Bull. Seism. Soc. Am. 72, 1843-1865.
- Boore, D. M., W. B. Joyner, and T. E. Fumal (1993). Estimation of response spectra
and peak accelerations from western North American earthquakes: An interim report,
U.S. Geol. Surv. Open-File Rept. 93-509, 72 pp.
- Brune, J. N. (1970). Tectonic stress and the spectra of seismic shear waves fromearthquakes, J. Geophys. Res., 75, no.26, 4997-5009.
- Campbell, K. W. (1981). Near-source attenuation of peak horizontal acceleration, Bull.
Seism. Soc. Am. 71, 2039-2070.
- Cocco, M. and J Boatwright (2004).The envelopes of acceleration time histories, Bull.
Seism. Soc. Am. 83, 1095-1114.
- Clinton, J. F. (2004). Modern digital seismology : instrumentation, and small amplitude studies in the engineering. Ph.D. thesis, Department of Civil Engineering,
California Institute of Technology.
- Cua, G. (2005). Creating the Virtual Seismologist: developments in earthquake early
warning and ground motion characterization. Ph.D. thesis, Department of Civil Engineering, California Institute of Technology.
- Cua, G. and T. H. Heaton (in press). The Virtual Seismologist method: a Bayesian
approach to seismic early warning. In ”Seismic early warning”, eds: P. Gasparini, J.
Zschau, Spring Heidelberg.

M. Yamada and T. Heaton

13- Hanks, T. C. and R. K. McGuire (1981). The Character of High-frequency Strong
Ground Motion, Bull. Seism. Soc. Am. 71, 2071-2095.
- Heaton, T. H. and S. H. Hartzell (1989). Estimation of Strong Ground Motions from
Hypothetical Earthquakes on the Cascadia Subduction Zone, Pacific Northwest, Pageoph 129, 131-201- Ji, C., D. V. Helmberger, D. J. Wald and K. F. Ma (2003). Slip history and dynamic implication of 1999 Chi-Chi earthquake. J. Geophys. Res., 108(B9), 2412,
doi:10.1029/2002JB001764.
- Kanamori, H. and P. C. Jennings (1978). Determination of local magnitude, ML, fromstrong motion accelerograms, Bull. Seism. Soc. Am. 68, 471-485.
- Lee, W. H. K., T. C. Shin, K. W. Kuo, K. C. Chen, and C. F. Wu (2001). CWB
Free-Field Strong-Motion Data from the 21 September Chi-Chi, Taiwan, Earthquake,
Bull. Seism. Soc. Am. 91, no.5, 1370-1376.
- Liu, H. L., and D. V. Helmberger (1985). The 23:19 aftershock of the 15 October 1979
Imperial Valley earthquake: more evidence for an asperity, Bull. Seism. Soc. Am.75,689-708- Ma, K. F., J. H. Wang, and D. P. Zhao (1996). Three-dimensional seismic velocity
structure of the crust and uppermost mantle beneath Taiwan, J. Phys. Earth, 44,
85-105.
- Nakamura, Y. (1988). On the urgent earthquake detection and alarm system (UrEDAS), Proc. of the 9th World Conference on Earthquake Engineering, Tokyo-Kyoto,
Japan. VII. 673-678.
- Nakamura, Y. (2001). Fault movement recorded by UrEDAS during the 2001 Geiyo
Earthquake, Japan, Proc. of the 26th JSCE Earthquake Engineering Symposium,
Sapporo, Japan. 1409-1412 (in Japanese).
- Odaka, T., K. Ashiya, S. Tsukada, S. Sato, K. Otake, and D. Nozaka (2003). A
New Method of Quickly Estimating Epicentral Distance and Magnitude from a Single
Seismic Record, Bull. Seism. Soc. Am. 93, no.1, 526-532- Sambridge, M. (1999a). Geophysical Inversion with a Neighbourhood Algorithm-I.
Searching a Parameter Space. Geophys. J. Int. 138, 479-494.
- Sambridge, M. (1999b). Geophysical Inversion with a Neighbourhood Algorithm-II.
Appraising the Ensemble. Geophys. J. Int. 138, 727-746.
- Shin, T. C. and T. L. Teng (2001). An Overview of the 1999 Chi-Chi, Taiwan, Earthquake, Bull. Seism. Soc. Am. 91, no.5, 895-913.
- Wessel, P., and W. H. F. Smith (1991). Free software helps map and display data,
EOS. Trans. AGU 72, 441, 445-446.
- Wu, Y. M. and H. Kanamori (2005). Experiment on an onsite early warning method
for the Taiwan Early Warning System, Bull. Seism. Soc. Am. 95, 347-353.
- Yamada, M. (2007). Early warning for earthquakes with large rupture dimension.
Ph.D. thesis, Department of Civil Engineering, California Institute of Technology.
- Yamada, M., T. H. Heaton and J. L. Beck (2007). Real-time Estimation of Fault
Rupture Extent Using Near-source versus Far-source Classification, Bull. Seism. Soc.
Am. (in press).

14

M. Yamada and T. Heaton

Figures and Tables
Thickness(km)
1.0
3.0
5.0
4.0
4.0
8.0
5.0
Half-space

Vp (km/s)
3.50
3.78
5.04
5.71
6.05
6.44
6.83
7.28

Vs (km/s)
2.00
2.20
3.03
3.26
3.47
3.72
3.99
4.21

Table 1: P-wave and S-wave velocity model in central Taiwan (Ma et al., 1996).

15

M. Yamada and T. Heaton

24˚45'
nf
au
ihta

0

50

Sh

24˚30'

lt

km
StationEpicenter
Subsource

7
6
5

24˚15'
4

Taichung

32

23˚45'

C024

Chelungpu fault

24˚00'

Hualian

1
0
−1
−2

−3
−4

23˚30'

23˚15'
120˚00'Chiayi120˚15'120˚30'120˚45'121˚00'121˚15'

121˚30'121˚45'122˚00'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 thispaper. The area within 50 km 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 theVoronoi cell for the station.

M. Yamada and T. Heaton

16

Figure 2: Observed envelope for accelerogram and P-wave and S-wave envelopes for the groundmotion model defined in equation (5) (Cua, 2005).

17

N1

=4

fau
lt

M. Yamada and T. Heaton

azimuthalangle
stationN2=3

epicenterstationsubsource

stationFigure 3: Schematic diagram of the multiple source model. The fault rupture is assumed topropagate from the epicenter at the constant velocity vR . The fault is parameterized by θ, N1, andN2, where θ is the azimuthal angle of the fault, N1 and N2 are the number of subsources north andsouth of the epicenter, respectively. The ground motion at a station is expressed as a combinationof the envelope from each subsource.
acceleration (cm/s2)

150
Predicted (VS−FS)
100

50

0
0

20

40

60

80

100

acceleration (cm/s2)

250
Observed
Predicted (VS−PS)
Predicted (VS−FS)

200
150
100
50
0
0

20

40

60

80

100

time (sec)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 magnitude 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 (insolid black line).

18

M. Yamada and T. Heaton

25˚00'
km
24˚45'

0

20

T147

40

60

Ilan

T026
T095

StationEpicenter
Subsource
Predicted
Observed

24˚30'

T081

T085
I063

T035
T045
T038

I064
I067

T128

H045

24˚15'

T068

T059

T088

Taichung

H056
H026

T111

T067
T072

24˚00'

H019
Hualian

T115
T119

T076

23˚45'

T084

H059

C024
C082
H033

C104

C093

C080

23˚30'

C034
C008

C074

Chiayi

H038

C041

C012
C107

23˚15'N001
C052

23˚00'
120˚00'

H041

K001

C099
C116

K054

N041
N014

C061

Tainan120˚15'120˚30'120˚45'121˚00'121˚15'

121˚30'121˚45'122˚00'60 secs.

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 fromthe best fit model are shown in a small yellow circles. The area within 50 km and 100 km from theepicenter are shown by large circles. Only characteristic records of the stations are shown in thisfigure.

19

M. Yamada and T. Heaton

25˚00'
km
24˚45'

0

20

T147

40

60

Ilan

T026
T095

StationEpicenter
Subsource
Predicted
Observed

24˚30'

T081

T085
I063

T035
T045
T038

I064
I067

T128

H045

24˚15'

T068

T059

T088

Taichung

H056
H026

T111

T067
T072

24˚00'

H019
Hualian

T115
T119

T076

23˚45'

T084

H059

C024
C082
H033

C104

C093

C080

23˚30'C034

C074

Chiayi

H038

C041

C012
C107

23˚15'N001
C052

23˚00'
120˚00'

H041

K001

C099
C116

K054

N041
N014

C061

Tainan120˚15'120˚30'120˚45'121˚00'121˚15'

121˚30'121˚45'122˚00'60 secs.

Time (s) from event
Figure 6: Predicted and observed envelopes in vertical component. The red and black lines are thepredicted and observed envelopes, respectively. The locations of the subsources estimated from thebest fit model are shown in a small yellow circles. The area within 50 km and 100 km from theepicenter are shown by large circles. Only characteristic records of the stations are shown in thisfigure.

20

M. Yamada and T. Heaton

azimuthal angle θ (deg)

100

50

0
full dataset (239 stations)subset A (126 stations)subset B (56 stations)

−50
0

10

20

30
time(s)

40

50

60

50

60

N1 (north)

Vr=2.0 km/s

10
50(south) N2

Number of the subsources

15

−5
−10
−15
0

full dataset (239 stations)subset A (126 stations)subset B (56 stations)

10

20

30
time(s)40

Figure 7: Time series of the estimated parameters, θ, N1, and N2 for each dataset. Time is relative
to the origin. The parameters are computed at each second using only the data available at thattime. 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 for N1 and N2. The solid
thin lines are the upper limits for N1 and N2 for the rupture velocity 2 km/s.

21

M. Yamada and T. Heaton

azimuthal angle θ (deg)

100

50

0

−50

−100
0

Vr = 2.0km/s
Vr = 2.5km/s
Vr = 3.0km/s
Vr = 3.5km/s

10

20

30
time(s)

40

50

60

50

60

N1 (north)Vr=3.5 km/s

10

Vr=2.0 km/s

50(south) N2

Number of the subsources

15

−5
−10

Vr = 2.0km/s
Vr = 2.5km/s
Vr = 3.0km/s
Vr = 3.5km/s

−15
0

10

20

30
time(s)40

Figure 8: The estimated parameters, N1, and N2 for different rupture velocities. The solid thinlines are the upper limits for N1 and N2 for the rupture velocity 2 km/s and 3.5 km/s. The brokenlines 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.

22

M. Yamada and T. Heaton

60

misfit

55
50

55

45

50

40
0

45
2
4
6
8
N1

−30 0
−60
10 −90θ (deg)

30

60

90

Figure 9: Error surface of θ and N1 at the fixed N2 = 4 at 60 seconds after the origin time. Since
the surface is peaked around θ = 0, it is easy to converge in θ. However, the optimal N1 will change
easily depending on the misfit function (see equation (3)).

65
misfit

60

60

55
50

55

45

50

40
0

45
2
4
6
8
N2

10 0

2

4

6

8

10

N1

Figure 10: Error surface of N1 and N2 at the fixed θ = 17 at 60 seconds after the origin time.
Since the surface is smooth in both N1 and N2 direction, the optimal solution is sensitive to a small
disturbance.

23

M. Yamada and T. Heaton

10

10

10

T = 10sec

T = 20sec

T = 30sec

6

6

6
N2

8

N2

8

N2

8

4

4

4

2

2

2

0

0

2

4

6

8

0

10

0

2

4

N1

6

8

0

10

10

4

4

4

2

2

2

N1

8

10

10

N2

6

N2

6

N2

6

6

8

T = 60sec
8

4

6

10

8

2

4

T = 50sec

8

0

2

N1

10
T = 40sec

0

0

N1

0

0

2

4

6
N1

8

10

0

0

2

4

6

8

10

N1

Figure 11: Contour maps of the error surface of N1 and N2 at the fixed θ = 17. The maps areshown in 10 second 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/s.