Bulletin of the Seismological Society of America, 91, 5, pp. 1112–1127, October 2001

The Importance of the Dynamic Source Effects on Strong Ground Motion
during the 1999 Chi-Chi, Taiwan, Earthquake: Brief Interpretationof the Damage Distribution on Buildings
by L. A. Dalguer, K. Irikura, J. D. Riera and H. C. Chiu

Abstract The 1999 Chi-Chi, Taiwan, earthquake, that originated on a low-angle
reverse fault, showed complexity and uncommon characteristics. The records showthat the hanging-wall side is characterized by larger particle motions than the footwall, and the ground motion is stronger in the northern part than in the southern partof the causative fault. Although the strongest ground motion occurred near the northern part of the trace, structural damage was heavier in the southern part. In order to
get a better understanding of the complex damage distribution caused by this earthquake, the dynamic rupture process was numerically simulated. Because of the differences between the observed features of the rupture process in the northern andsouthern parts of the fault, each part was modeled independently by using a 2D
discrete element model (DEM). The principal results of the simulation show that thevelocity ground motions in the northern part, in the frequency range of 0.5–2 Hz(natural frequency range of standard structures), are small near the surface break,thus, light structural damage might be predicted near the surface rupture. Moreover,in the northern part the fault rupture propagation reaches the surface with a very
slow velocity (about 1.2 km/sec); however, in the southern part the rupture propagation reaches the surface with higher velocity (about 3.0 km/sec). These differences
between the models could explain why the ground motion near the surface rupturein the northern part caused less damage in structures than the ground motion in thesouthern part.
Introduction
Due to its complexity and uncommon characteristics,
the Chi-Chi, Taiwan, earthquake (Mw 7.6) of 20 September
1999 had a great impact on the international community ofscientists and engineers devoted to seismology and earthquake engineering. The earthquake originated on a lowangle reverse fault with a strike of nearly N5⬚E and a dip
between 25⬚ and 36⬚ (Shin et al., 2000). The rupture of thecausative fault reached the surface and propagated along
about 80 km, starting in the south and extending northward
on the Chelongpu fault, as shown in Figure 1. Spectacularhorizontal displacements up to 9.0 m and vertical offsets of
1.0 to 8.0 m were registered along the surface rupture. The
earthquake was one of the most destructive and largest earthquakes of this century in Taiwan. The death toll was 2333,
with 10002 people injured (Shin et al., 2000). The epicenterof this shallow earthquake was located near the small townof Chi-Chi (23.87N, 120.75E), 150 km (90 miles) south of
Taipei. In order to analyze the damage distribution causedby the earthquake, the northern and southern zones are defined as shown in Figure 2a and b. Although the stations

records show that the strongest ground motion and the largest displacements were located in the northern part of the
Chelongpu fault, structural damage in that area appears tobe less than in the southern portion. In fact, only buildings
directly on the surface rupture had severe damage in thenorthern zone. In general, the most severe damage on buildings caused by the ground shaking was registered in the
hanging wall and was greater in the southern part than in thenorthern portion (Fig. 2a). From the distribution of damaged
structures, shown in Figure 2a, and the population distribution along the causative fault, shown in Figure 2b, we canobserve that the northern part of the fault (except for its final
portion) has higher population density than the southern part(near the epicenter area); however, damage is larger in thesouthern part. Such difference in the damage distribution canbe inferred from Figure 3, which shows a comparison of thespectral pseudovelocity from records of stations located in
the north (TCU052) and south (TCU129 and TCU089). The
peak velocity for frequencies less than 1.0 Hz is larger in
the northern part. On the contrary, for frequencies greater

1112

The Importance of the Dynamic Source Effects on Strong Ground Motion during the 1999 Chi-Chi, Taiwan, Earthquake

1113

ground-motion records. We showed that this kind of fault
generates large differences between the near-source groundmotions on the hanging wall and on the footwall. The groundmotions on the hanging wall are larger than in the footwall,and the numerical simulation suggests that such a difference
is principally caused by the asymmetric geometry of the hanging wall and footwall. For this earthquake, where the ruptureof the fault reaches the surface, the effect of the asymmetryon the ground motion is considerable. The characteristics ofthis kind of earthquakes (dipping fault) were also shown by
Mikumo and Miyatake (1993) in their investigation of the
dynamic rupture process of the 1961 Kita-Mino earthquake
in central Japan and by previous theoretical dynamic simulation of dipping fault, for example Nielsen (1998), Shi etal. (1998), Oglesby et al. (1998, 2000), as well as that observed in the foam rubber experiment of a thrust fault presented by Brune (1996).

Dynamic Model

Figure 1.

Location of the surface rupture of the
Chelongpu fault; stations records used for comparisonand sections of the northern and southern model.

than 1.0 Hz, the peak velocity ground motion is larger in thesouthern part. This means that the northern region generated
stronger ground motion at low frequencies than did thesouthern region. Figure 2a also shows the highest concentration of damage in a small area at the end of the fault in
the northern part, where the rupture of the fault stopped,
turning suddenly to the east. This region of the fault is outof our scope because of its particular complexity, escaping
from the purpose of this article to explain the general damage
distribution caused by the earthquake.
In order to better understand the complex damage distribution caused by this earthquake, we performed a numerical simulation of the dynamic rupture process. Since thefaulting process appears to have been nearly pure thrustalong the various fault segments, a 2D discrete element
model (DEM) was employed to perform a dynamic simulation of the rupture process and of the near-fault ground motion. On account of the differences in the observed featuresof the rupture process in both parts of the fault, each partwas modeled independently. The dynamic rupture processsimulation of the southern part was already presented by theauthors in a previous article (Dalguer et al., 2001). In thisarticle we simulated the mechanism of the reverse fault near
the hypocentral area, as well as the near-source ground motion, the results obtained closely matching the observed

Several attempts to perform numerical analyses of the
dynamic rupture processes of a fault have been described in
the literature. The pioneering work of Kostrov (1966) simulated the spontaneous propagation of an anti-plane shear
crack. Das and Aki (1977) and Andrews (1976), using theslip-weakening model as a friction law of the fault, simulated
spontaneous rupture propagation of an in-plane shear crack.
Subsequently, the rupture process of the fault was simulatedwith more sophisticated models (e.g., Mikumo and Miyatake, 1978; Day, 1982a,b; Virieux and Madariaga, 1982;
Cochard and Madariaga, 1994, Fukuyama and Madariaga,
1998; Inoue and Miyatake, 1998; Madariaga et al., 1998).
Dynamic models are frequently used to study the physics ofearthquakes, as related to the rupture process of the fault.
Although few efforts were devoted to the simulation ofground motion based on dynamic models, some recent contributions should be cited. Olsen et al. (1997) simulated therupture process and near field ground motion of the 1992
Landers (California) earthquake using a finite-difference
method in the frequency range 0.1 to 0.5 Hz. Inoue and
Miyatake (1998) simulated theoretical strong ground motion
generated from the rupture process on a shallow strike-slip
fault using a 3D finite-difference method in frequencies up
to 2.0 Hz.
In this article, we simulate the rupture dynamic of apreexisting fault and near-source ground motion, solving theelastodynamic equation of motion using the discrete element
method (DEM). The DEM is widely employed in engineering
to designate lumped mass models in a truss arrangement, as
opposed to finite element method (FEM) models that may
also consist of lumped masses, but normally require to
mount a full stiffness matrix for response determination. The
term has also been used for models of solids consisting ofassemblies of discrete elements, such as spheres in elastic
contact, employed in the analysis of perforation or penetration of concrete or rock. It should be noted that the desig-

1114

L. A. Dalguer, K. Irikura, J. D. Riera, and H. C. Chiu

Figure 2.

(a) Damage distribution of buildings in Nantou and Taichung area: orange
dots show completely damaged areas and green dots show severely damaged areas(Architecture & Building Research Institute, Ministry of Interior, ROC, Taipei, Taiwan); (b) Population distribution along the causative fault (Tsai and Huang, 2000). The
red line represent the Chelungpu fault (causative fault of the 1999 Chi-Chi earthquake),and the purple line in part a is the Shuangtung fault.

Figure 3.

Comparison of pseudo–spectral velocity for records of stations at thenorthern and southern parts.

nation lattice models, common in physics, may be more adequate, although it omits reference to a fundamental propertyof the approach, which is the lumped-mass representation.
In the present DEM formulation, orthotropic solids are represented by a three-dimensional periodic trusslike structureusing cubic elements as shown in Figure 4. This model isbased on earlier developments in aeronautical engineeringin which, for purposes of structural analysis, it is often necessary to establish the equivalence between trusslike structural systems and a continuous medium. Nayfeh and Hefsy(1978) established the equivalence requirements between the
cubic arrangement shown in Figure 4 and an orthotropic
elastic medium. The method leads to results that converge
to solutions for a linear elastic continuum in dynamic problems. Riera and Rocha (1991) used the approach in fracturestudies, Doz and Riera (1995) employed the method tomodel the stick-and-slip motion along friction surfaces, Dalguer et al. (1999) evaluated the foreshock and periodicity ofearthquakes, and Dalguer et al. (2001) modeled a fault dynamic rupture using the slip-weakening friction model toThe Importance of the Dynamic Source Effects on Strong Ground Motion during the 1999 Chi-Chi, Taiwan, Earthquake

1115

that reflects the random distribution of the orientation of constituent elements. Details of the calculation of the equivalentcross-sectional axial stiffness of the normal (AEn) and diagonal (AEd) elements for a cubic lattice array given byequation (1) and (2) respectively are presented in the Appendix.
In the discrete dynamic model, masses are concentrated
at nodal points. As shown in Figure 4a, solids are represented
as an array of normal and diagonal bars linking lumped
nodal masses. The dynamic analysis is performed using explicit numerical integration in the time domain. At each stepof integration a nodal equilibrium represented by equation(3) is solved by the central finite-differences scheme.
müi Ⳮ cu̇i ⳱ fi ,Figure 4.

Numerical model used for the dynamic
simulation (DEM): (a) basic cubic module, (b) generation of prismatic body, and (c) representation of aplate-strain state (no z displacements).

simulate the rupture process of the hypocentral area of the1999 Chi-Chi (Taiwan) earthquake. Also Mora and Place(1994), Shi et al. (1998), and Rimal (1992) used a 2D latticemodel similar to the DEM to simulate dynamic rupture ofearthquake faulting.
In the case of an isotropic elastic material, the crosssectional axial stiffness of the longitudinal bars in the equivalent discrete model is given by:
AEn ⳱ ␾EL2, bar length ⳱ L,

(1)

while for the diagonal bars,
AEd ⳱ 2d␾EL2冪3

, bar length ⳱ 冪3

L
,
2

(2)

where for approximately isotropic solids, namely, solids
with equal stiffness in the three orthogonal directions, ␾ ⳱(9 Ⳮ 8 d)/(18 Ⳮ 24 d), d ⳱ 9v(4 ⳮ 8v), v is the Poisson’s
ratio, and E is the Young’s modules of the material. For othersituations, for example layered (orthotropic) rocks, the aforementioned constants take on other values (See Nayfeh andHefsy, 1978). It should be stressed that no lattice or trusslike
model can exactly represent a locally isotropic continuum,
and for that matter it can also be argued that no locally isotropic continuum exists. Isotropy in solids is a bulk property(3)

where m denotes the nodal mass, c is the damping constant,ui , u̇i , üi denote a component displacement, velocity, andacceleration, respectively, of the nodal coordinates vector,
and fi denotes a component of the resultant forces at one
nodal point including elastic, external, and frictional forcesin direction i of the motion. In the current model, only thenodal points that coincide with the preexisting fault, once it
breaks, are under frictional force governed by any predefined
friction law. The damping constant c was assumed to be
proportional to the rigidity k of the bars of every cubic element, that is c ⳱ df k, where df was assumed to be 0.005. It
is approximately a critical damping ratio (v) less or equal to0.045.
Equation (3) represents the equation of motion of a discrete point in the continuum. When this point is on the preexisting fault, the fault parallel component of the resultant
force (fi) is governed by the constitutive relation on the fault.

Friction Law on the Fault
Laboratory experiments on rock (e.g., Dieterich, 1979;
Ohnaka et al., 1987, Ruina, 1983) lead to slip- and/or ratedependent friction models. In the present model we adopt
the simple slip-weakening friction model that was first proposed by Ida (1972) and extensively used for dynamic simulation of fault rupture processes (e.g., Andrews, 1976; Day,
1982b; Olsen et al., 1997; Fukuyama and Madariaga, 1998;
Harris and Day, 1999). The slip-weakening friction modelis schematically represented in Figure 5. The shear force
calculated from the resultant force, fi, of equation (3) could
be expressed by the shear stress, s. The following is the relation between the shear stress, s, and the slip of the fault, u.
s ⬍ su ,for u ⳱ 0,su ⳮ sfs⳱ⳮu Ⳮ su for 0 ⬍ u ⬍ Dc ,Dcs ⳱ sffor u ⱖ Dc ,(4)

where su is the critical stress or the upper yield point, sf isthe final stress or the residual stress, which is considered as

1116

L. A. Dalguer, K. Irikura, J. D. Riera, and H. C. Chiu

Figure 5.

The slip-weakening friction model.

the dynamic friction stress level, and Dc is the slip required
for stress to drop to its dynamic friction level. We assume
that there is not back slip on the fault, which means that theslip velocity is always greater or equal to zero.
In order to verify the adequacy of our approach, using
the DEM, to simulate a dynamic rupture process, we analyzed the spontaneous inplane rupture with the slip weakening law employed by Andrews (1976). In this problem,
plane strain is assumed. The crack plane is the x-y plane, inwhich the crack propagates bilaterally in the x direction andextends indefinitely in the y direction. The medium is thus
an infinite, homogeneous, isotropic, and linearly elastic
crack plane. As indicated before, the rupture assumes a slipweakening law. When the crack propagates, it will not stop.
For the numerical computations, the parameters are normalized as follow (Andrews, 1976): Shear stress along thecrack plane s⬘ ⳱ s/Ds, x axis parallel to the crack planex⬘ ⳱ x/Lc, time t⬘ ⳱ tb/Lc; slip: u⬘ ⳱ ul/LcDs; slip velocity:v⬘ ⳱ vl/bDs, where b is the S-wave velocity, l is the shear
rigidity, Ds is the stress drop, and Lc is the critical half-lengthof a Griffith crack in plane strain derived by Andrews(1976),Lc ⳱

8l(k Ⳮ l)G,p(k Ⳮ 2l)(so ⳮ sf)2(5)

where k and l are the Lamé constants, and G is the effective
fracture surface energy given byG⳱1(s ⳮ sf)Dc .
4 u

(6)

We assume that that Poisson’s coefficient is 0.25, so
␣/b ⳱ 冪3, where ␣ is the P-wave velocity. These nondimensional quantities are equivalent to assuming that l ⳱ 1,
Ds ⳱ 1, b ⳱ 1, ␣ ⳱ 冪3, Lc ⳱ 1, density q ⳱ 1. The
calculations were performed with a grid size Dx ⳱ 0.1Lc(length of the side of one cubic element) and (su ⳮ so)/Dr⳱ 0.8. Using equations (5) and (6), Dc ⳱ 1.31.
The numerical solution is obtained for the near-fieldelastodynamic motion coupled to frictional sliding on a prescribed crack plane. Initially the stress distribution along thefault is in the initial stress level (so) shown in Figure 5, therupture is initiated artificially by imposing a stress drop to
propagates at least as fast as 0.5b, which leads to initialstresses along the fault that increase monotonically without
any relative slipping along the fault. Eventually, the interface
shear stress (s) at a point exceeds the local shear strength(critical stress level su) and slip at a node occurs, governedby the slip weakening model shown in Figure 5 and represented by equation (4). Considering that the seismic radiation depends only on the stress change (stress drop) during
the earthquake, and not on the total stress, we assume thatthe initial stress (so) along the fault is at its zero level. Therefore, the necessary parameters required to simulate the rupture process governed by the slip-weakening friction model
are the strength excess, stress drop, and critical slip.
Figures 6 shows the results of our numerical simulation
for a theoretical spontaneous inplane rupture problem analyzed by Andrews (1976). The results are very consistent
with the solution presented by Andrews (1976). Figure 6ashows the space–time distribution of rupture. The region between the two solid lines is the rupture front, where slip
velocity is nonzero and stress drop is incomplete. Figure 6bshows the slip velocity as a function of position at the same
instant, at dimensionless time bt/Lc ⳱ 12.36. Figure 6c–fshows the shear stress and slip as a function of position on
the crack plane at dimensionless time bt/Lc ⳱ 8.07, 10.38,
12.36, and 14.34, respectively. The peak of the shear stresson the rupture front as well as the secondary peak associatedwith the S waves are very precisely described by the DEM.

Assumptions in the Numerical Simulation to
Simulate the 1999 Chi-Chi Earthquake
Because of the differences in the observed features of
the rupture process in the northern and southern parts of thecausative fault of the 1999 Chi-Chi earthquake, each partwas modeled independently. The problem is tackled in aplane-strain condition. The location of the two models sections along the surface rupture is shown in Figure 1. The
first model (southern part) is near the epicentral area, andthe second model (northern part) is near the TCU052 station.
The parameters used for the dynamic simulation and the geometry of the two fault models are shown in Figure 7. The
fault model of the southern part is the same used by Dalguer
et al. (2001). Both models share the following common assumptions:
1. There is a surface sedimentary layer with a depth of 4 km
characterized by a set of P-wave velocity (4.3 km/sec),
S-wave velocity (2.5 km/sec), density 2500 kg/m3, corresponding to Young’s modulus 3.9 ⳯ 1010 N/m2, shear
modulus 1.56 ⳯ 1010 N/m2 and Poisson’s ratio 0.25. The
basement (seismogenic zone) is a homogeneous medium
with P-wave velocity (6.1 km/sec), S-wave velocity (3.5km/sec), density 2700 kg/m3, corresponding to Young’s

The Importance of the Dynamic Source Effects on Strong Ground Motion during the 1999 Chi-Chi, Taiwan, Earthquake

1117

Figure 6. Dynamic solution, resolved by the DEM, of the spontaneous inplane rupture problem presented by Andrews (1976). (a) Space–time of rupture propagation.
Region between the two solid lines is the rupture front, where slip velocity is nonzero
and stress drop is incomplete. Dashed line labeled P, S, and R represents the wave frontof the compressional, shear, and Rayleigh waves, respectively, as a reference (Rayleighwave velocity ⳱ 0.9194b). (b) Dimensionless slip velocity, v⬘ ⳱ vl/bDs, as a functionof position on the crack at time the dimensionless time bt/Lc ⳱ 8.07. Shear stress andslip as a function of position on the crack plane at dimensionless time bt/Lc ⳱ (c) 8.07,(d) 10.38, (e) 12.36, and (f) 14.34, respectively. Heavier solid curve is dimensionlessslip function divided by 10, lu/10Lc (Ds); lighter solid line is the dimensionless changeof shear stress, (s/Ds).

modulus 8.37 ⳯ 1010 N/m2, shear modulus 3.35 ⳯ 1010
N/m2, and Poisson’s ratio 0.25.
2. The preliminary report of the 1999 Chi-Chi earthquake(Shin et al., 2000), suggests a focal depth of 7–11 km
and dip of 25⬚–36⬚. We selected a dip of 33⬚41⬘ in orderto have an integer element size of the DEM to model thefree surface. The hypocenter is located at a depth of
8.5 km.
3. The slip-weakening friction model is adopted as the constitutive relation for the fault.
4. The stress drop along the fault plane in the shallow surface layer is negligible.
5. The ultimate stress, that is, the strength excess on the fault
surface in the shallow surface layer increases linearlywith depth and, in order to avoid any fault opening, weapplied a normal stress along the fault, equivalent to thestrength excess.
The choice of the parameters used for the simulationof the dynamic rupture of a fault is a delicate issue, still
subject to debate. The author’s experience, nevertheless,suggests that, in general, it is very unlikely that essentiallycorrect simulated response will result from essentially incorrect models or parameters. For the assumption of a slipweakening friction law, we need to define the stress drop,
the strength excess, and the critical slip, Dc, along the fault.
The stress drop can roughly be defined from the results of
the kinematic inversion of ground motion, so we can guessapproximate asperities area (we assume that asperity is azone of higher stress drop than surrounding areas). The

1118

L. A. Dalguer, K. Irikura, J. D. Riera, and H. C. Chiu

Figure 7.

Fault models and parameters distribution used for the dynamic simulation:(a) southern model, (b) northern model.

choice of Dc and the strength is quite difficult, since there is
no way to get them from observations, and there may bevarious combinations that fit the data. This is in accordancewith the conclusion from Guatteri and Spudich (2000), who
indicate that there is no unique solution to select Dc andstrength excess. But in general, we can get some information
on Dc from several laboratory experiments and field analyses. Ohnaka (1992) reported that Dc is constant within thebrittle layer, Marone and Kilgore (1993) proposed that Dc islarger when a fault gouge zone is thicker, so, the rupture in
the shallow crust is stable due to the presence of a large
amount of gauge materials between the fault zone. And also,
from the report of Ide and Takeo (1997) who study the stresschange on the fault during the dynamic rupture propagationin the 1995 Kobe earthquake through seismic inversion, Dc
is more than 1.0 m in the shallow crust and about 10 cm inthe brittle layer. From our numerical experiments, we observed that the frequency content of the ground motion isrelated to the rupture propagation velocity, a high rupture
velocity producing higher frequencies than a slow velocity.
The rupture velocity is related to Dc and the strength excess,that is, small Dc and strength excess produce higher rupture
velocity than large Dc and strength excess. In view of thepreceding arguments, the modeling procedure in the article
was as follows: from kinematical inversion results, we canroughly get the asperities areas and the stress drops. The
preliminary results of kinematical inversions (e.g., Iwata etal., 2000; Ma and Mori, 2000; Yagi and Kikuchi, 2000)
show the existence of larger asperity located in the northernpart and relatively two or three small asperities in the southern part. The Dc and strength excess were selected in order
to fit the damage pattern, that is, define combination of Dc
and strength excess in order to explain in a general way thedifference of the ground motion between the northern andsouthern zone of the fault. Considering these observation,and the arguments given previously, we start from the assumption that Dc is larger near the free surface than at higher
depths and that it is larger in the northern region that in thesouthern region. Therefore, in the model for the southernpart, the existence of three asperities with relatively small
widths (about 7 km) and stress drop between 1.5 MPa and3 MPa were used in the basement underlying the sediments,
and the critical slip 0.5 m in the surface layer and 0.1 m inthe deep part (Fig. 7a). The northern part, on the other hand,is assumed to have just one asperity with larger width (15km) and higher stress-drop (between 1.5 MPa and 8 MPa),
the critical slip is 2.5 m in the surface layer and of 2.0 m inthe deep part (Figure 7b).
The fault models are constructed taking into account thetwo sides of the fault. We assume that the fault has a thickness equivalent to the size of one cubic cell. Once the fault
breaks, the linkage between the two surfaces of the fault isbroken, being governed by the friction law. The shear
stresses on both sides of the fault are equal in magnitude but
act in opposite directions. The models used for the simulation (Fig. 7) include a preexisting fault 40 km wide in a solidof 66 km ⳯ 28.25 km. 29832 cubic cells with 0.25 km longsides were needed to construct the 2D model space bounded
by a free surface and three artificial boundaries that simulate
the semiinfinite space. Viscous dampers were introduced on
the artificial boundaries in order to avoid the reflection ofwaves from the edges of the domain of computation. The
boundary conditions at the surface were specified as zero
shear and normal stresses. During the rupture process, only
the inertial stresses caused by the dynamic motion act on thesurface. For the numerical integration of the equation of motion, a time step of 0.05 sec was used.

The Importance of the Dynamic Source Effects on Strong Ground Motion during the 1999 Chi-Chi, Taiwan, Earthquake

Simulation Results and Discussion
The two fault models present similar characteristics, thatis, the results of the dynamic rupture process of the twomodels illustrate the dynamic effects of a thrust faulting on
the ground motion. The main characteristics of dipping
faults that breaks the surface are that the slip of the fault andthe ground motion is larger in the hanging wall than in the
footwall. As we analyzed in a previous article, where we
simulated the dynamic rupture process of the southern part
(Dalguer et al., 2001), the large difference in the near-source
ground motion between the hanging wall and footwall iscaused mainly by the asymmetric geometry of the hanging
wall and footwall. We also concluded that the effect on the
ground motion of the rupture reaching the surface is strong
near the surface break of the fault, because the radiatedwaves are trapped in the hanging wall, resulting in the amplification of the ground motion in the wedge of the hanging
wall.
The slip and slip velocity distribution along the fault
shown in Figure 8a (southern model) and b (northern model)suggest that the total slip duration along the fault and thefinal slip are larger in the northern model than in the southern

Figure 8.

1119

model. The peak slip velocity in the asperities of the faultof the southern model is reached suddenly as a damped step
function. This is similar to the Kostrov’s pulse (Kostrov,1964) that was theoretically obtained in a self-similar shearcrack propagation. Nevertheless, although the northernmodel also shows this sudden increase of the slip velocity
in the beginning of the rupture, the peak value is reachedsmoothly after passing this sudden step. These results suggest that the rupture process in the southern model is rougher
than in the northern model, despite the fact that the slip velocity near the surface is larger in the northern model thanin the southern model, as shown in Figure 9, which shows
the space–time distribution of slip velocity of the southernand northern model.
The comparison of the rupture propagation (Figure 10)shows that the rupture of the northern model crosses theasperity with supersonic velocity and suddenly changes tosubsonic velocity when it leaves the asperity. The rupture
velocity toward the free surface reaches small values between 0.9 km/sec and 1.8 km/sec, being the smallest near
the free surface. Toward the deeper part, the rupture reaches
velocities equivalent to the S-wave velocity, after leaving

Slip and slip velocity along the fault: (a) southern model, (b) northern model.

1120

L. A. Dalguer, K. Irikura, J. D. Riera, and H. C. Chiu

Figure 9.

Space-time of slip velocity: (a) southern model, (b) northern model.

Figure 10. Comparison between the northern andsouthern model of the space–time of rupture propagation.

the asperity. On the other hand, the rupture of the southernmodel reaches velocity equivalent to the S-wave velocity
toward the free surface, but near the free surface the rupture
accelerates reaching supersonic velocity. Similarly, toward
the deeper zone, the rupture starts with subsonic velocity andincreases monotonically reaching supersonic velocity in the
deepest part of the fault. These results suggest that the rupture velocity in the northern model propagates slower thanthat of the southern model. The major difference is near thefree surface where the rupture velocity of the northern model
reaches in average 1.2 km/sec, while the southern model
reaches on average 3.0 km/sec.
The comparison of the ground motion on the surface

between the northern and southern model (in a frequency
range up to 2.00 Hz) is shown in Figures 11 and 12. Figure
11 shows the final displacements, peak velocity, and peakvelocity in a frequency range of 0.5 to 2.0 Hz. Figure 12
shows the waveforms of the velocity ground motion in afrequency range of 0.5 to 2.0 Hz. Figure 11a and h showsthat the final displacement and the peak velocity are larger
in the northern model than in the southern model. These
results show that the southern model predicts a vertical displacement around 2.0 m and a horizontal displacement
around 3.5 m on the hanging wall; on the other hand, thenorthern model predicts a vertical displacement around 4.0
m and a horizontal displacement around 6.5 m on the hanging wall. However, when the ground-motion velocities arefiltered in a frequency range of 0.5 to 2.0 Hz, as shown in
Figures 11c and 12, the peak velocities are larger in thesouthern model than in the northern model (Fig. 11c), andalso, the waveforms of Figure 12 show that the southernmodel contains shorter wavelengths than the northern model,
especially on the hanging wall. These results suggest thatthe northern model predicts stronger ground motions than
the southern model for lower frequency. On the other hand,in a higher frequency range, between 0.5 and 2.0 Hz (natural
frequency range of standard structures: 1–4 story), theground motion is stronger in the southern model. From theseresults we can conclude that although the northern model
presents stronger ground motions, the most severe damage
in structures should occur in the southern model, which ismore likely to severely excite the lower modes of standard
structures.
This differences of ground-motion intensity and predominant frequency between the northern and southernmodel may be due to the effects of the dynamic parametersof the source used for the simulation, that is, the resistanceof the fault (strength excess) and critical slip that are larger

The Importance of the Dynamic Source Effects on Strong Ground Motion during the 1999 Chi-Chi, Taiwan, Earthquake

1121

Figure 12. Comparison between the northern andsouthern model results of the waveforms of the velocity ground motion in a frequency range of 0.5–2.0
Hz.

Figure 11. Comparison of the final displacement
and peak velocity on the surface between the northernmodel and southern model: (a) final displacement, (b)peak velocity, (c) peak velocity in a frequency rangeof 0.5–2.0 Hz.

in the northern model than in the southern model. Whenthese parameters increase in magnitude, the rupture processof the fault is delayed, causing slower rupture velocity as
shown in Figure 10, where the rupture propagation of thenorthern model is slower than that on the southern model.
The critical slip also affects the beginning of the slipping
process. In this zone of the slip, called the breakdown zone,
the slip velocity reaches its peak as a damped step function
similar to the Kostrov’s pulse, and the magnitude of thispeak increases when the critical slip decreases. This effectcould be observed in Figure 8a (southern model) and b(northern model). These figures show that this peak is larger
in the southern model than in the northern model, eventhough the northern model fault has a larger stress drop. Thisis because the southern model has smaller critical slip than
the northern model. The effects of these parameters makethe southern model rougher in high frequencies than thenorthern model.
In order to observe the effects of the faulting when therupture breaks the surface, in Figure 13 we compare the final
displacements and the peak velocities on the surface fromtwo dynamic source models: one is the southern model
shown in Figure 7a in which the rupture reaches the ground1122

L. A. Dalguer, K. Irikura, J. D. Riera, and H. C. Chiu

Figure 13. Comparison between the southern model used to simulate the 1999 ChiChi earthquake in which the rupture breaks the surface (solid line), and another modelin which the rupture does not reach the surface (dashed line): (a) vertical and horizontal
component of the final displacement along the surface, and (b) vertical and horizontal
component of peak velocity along the surface.

surface (solid line), and the other one is a similar model inwhich the rupture is forced to stop 3 km before it reaches
the surface (dashed line). It can be clearly observed thatwhen the rupture reaches the surface there is a strong effecton the ground motions very close to the surface trace. The
ground motion is enhanced when the rupture breaks theground surface.

Figure 14. Comparison between the numericalsimulation and the observations (GPS data and stations records) of the final displacements along the surface near the epicentral area (southern model).

In order to validate the dynamic model studied here, we
compared our results of the southern model with the observations. Figure 14 shows that the southern model predicts
final vertical displacement of about 2.0 m and horizontal
displacement of about 3.3 m in the hanging wall that agree
satisfactorily with those obtained by the Global Positioning
System (GPS) data and the final displacements of the station
records. In Figure 15a–e we also compare the waveform of
the displacement and velocity ground motions of east–west
and vertical components recorded at five stations near thesurface rupture of the epicentral area (stations TCU084 andTCU089 on the hanging-wall side and TCU129, TCU116,
and TCU122 on the footwall side). We find that the main
characteristics of the recorded ground motion are adequately
reproduced, as observed in the simulated and observed time
histories. In the frequency range from 0.5 to 1.0 Hz, althoughobtained with a 2D model, the ground-motion simulation
qualitatively matches the observations. The simulations forstations TCU084 and TCU089, located on the hanging wall,

Figure 15. Comparison of the numerical simulation with stations records of horizontal (E–W) andvertical component of displacement, velocity, andvelocity in frequency range of 0.5–1.0 Hz: (a) for
TCU084 station, (b) for TCU089 station, (c) for
TCU129 station, (d) for TCU122 station, (e) for
TCU116 station, (f) TCU052 station. Dh and Dv are
the horizontal (east-west) and vertical components,
respectively, and Vh and Vv are the horizontal (east–west) and vertical component of velocity ground motion, respectively.

The Importance of the Dynamic Source Effects on Strong Ground Motion during the 1999 Chi-Chi, Taiwan, Earthquake

1123

1124

L. A. Dalguer, K. Irikura, J. D. Riera, and H. C. Chiu

fit the observations well, except for the horizontal component of the TCU084 station, which, according to Chiu
(2000), has noticeable site effects in the horizontal component. On the other hand, the simulations for stations
TCU129, TCU122, and TCU116, located on the footwallside, agree well with the observed motion only in the initialpart of the complete time histories. The ground motion on
the footwall side was maybe influenced by softer surfacelayer effects that are not taken into account in our simulation.
The simulation of the ground-motion records by meansof a 2D model is more feasible for the southern model because herein the 3D effects would be smaller in relation to
the geometry from the hypocenter to the site. Consequently
we can in such case perform a robust simulation using asimplified 2D model. On the other hand, we believe that thenorthern part was strongly affected by the 3D effects, rendering the 2D ground-motion simulation records highly
questionable. In this context, the dynamic rupture simulationof the northern model fault was carried out just to allow
some comparisons with the southern part, in connection withthe damage distribution, which is the objective of this article.
But even with the limitation of the northern model mentioned previously, we succeeded in simulating the groundmotion recorded at TCU052 station, as shown in Figure 15f,
in which the simulation closely fits the recorded motion. Thisis because this station is located very close to the fault on
the hanging wall, as shown in Figures 1 and 7b, where thedislocation of the fault is clearly recorded in the data.

Conclusions
Certainly there will be widespread disagreement on thepossibility of representing a complex earthquake such as the1999 Chi-Chi, Taiwan, earthquake by means of a simple 2Dmodel. We want to stress that our intention was not to fully
simulate the rupture process of the fault nor the resulting
ground motions, for which purpose a complete 3D model
would be needed, but to provide some interpretation of the
distribution of damage on structures caused by the earthquake, based on the dynamic rupture process.
Using the simplified 2D DEM, it was possible to show
some aspects of the dynamic source effects on the strong
ground motion of the 1999 Chi-Chi earthquake. The groundmotion simulation of the southern model has a satisfactory
agreement with the observed records. Also, even with thelimitation of the northern model to simulate ground-motion
records, the comparison of the simulation with the TCU052
station record is very good. These comparisons with the observations provide the validity of the results presented in thisarticle.
The results show that the velocity ground motions in thenorthern part (hanging wall) in the frequency range of 0.5–
2 Hz (natural frequency range of standard structures) aresmall near the surface break, thus, light structural damage
might be predicted near the surface rupture. The model ofthe southern part presents smaller displacements and smaller

velocity ground motions than the model at the northern part.
As discussed before, in spite of this global difference, theground motion in the southern part, in the frequency rangeof 0.5–2 Hz, is larger than in the north. Moreover, the fault
rupture propagation reaches to surface with a very slow velocity (about 1.2 km/sec) in the northern part. However, inthe southern part, it reaches the surface with considerable
higher velocity (a supersonic velocity of about 3.0km/sec).
The supersonic rupture velocity is a controversial issue. In
our simulation, supersonic velocity is observed near the freesurface in the southern model because we use small strength
excess and small Dc. These values were adopted in order togenerate higher frequency ground motion than in the northern model. In this context, the near source high frequency
ground motion actually observed in the strong-motion data
could be caused by possible supersonic rupture velocity inthe southern zone of the Chi-Chi earthquake.
These results suggest that the ground motion near thefault of the 1999 Chi-Chi earthquake was affected mainly by
dynamic source effects. These source effects predict stronger
ground motions in the north for lower frequencies. However,for higher frequencies, between 0.5 and 2.0 Hz, the groundmotion was predicted to be stronger in the southern model.
From these results we can conclude that even though thenorthern model shows stronger ground motion, the most severe damage to structures should be expected in the south
because it has a greater potential to severely excite the fundamental mode of most standard, low height structures. The
model presented in this article clearly shows that dynamic
effects of the source mechanism on strong ground motion
prediction are of fundamental importance for the assessmentof seismic hazard.

Acknowledgments
We are thankful to the Seismology Center, Central Weather Bureau(CWB), Taipei, Taiwan, for providing the Digital waveform files corresponding to the accelerograms of the stations used, and to Prof. Kuo-Liang
Wen from the Office of the National Science and Technology Program for
Hazards Mitigation, Taipei, Taiwan, for his cooperation in releasing information on structural damage and population distributions population distribution along the causative fault. We express our sincere thanks to thereviewers for constructive criticism and valuable suggestions. This research
was supported by a Grant-in-Aid for Science Research, No. 11792026, from
the Ministry of Education, Science, Sports and Culture of Japan, and by
the National Council for Scientific and Technological Development (CNPq)
of Brazil.

References
Andrews, D. J. (1976). Rupture velocity of plane-strain shear cracks, J.
Geophys. Res. 81, 5679–5687.
Brune, J. N. (1996). Particle motion in a physical model of shallow angle
thrust faulting, Proc. Indian. Cad. Sci. 105, 197–206.
Chiu, Hung-Chie (2000). Strong ground motion records from the 1999 ChiChi earthquake, unpublished document, Institute of Earth Science,
Academia Sinica, Taipei, Taiwan.
Cochard, A., and R. Madariaga (1994). Dynamic Faulting Under RateDependent Friction, PAGEOPH 142, 419–445.

The Importance of the Dynamic Source Effects on Strong Ground Motion during the 1999 Chi-Chi, Taiwan, Earthquake

1125

Dalguer, L. A., J. D. Riera, and K. Irikura (1999). Simulation of seismic
excitation using a stick-slip source mechanism, Transaction, 15th International Conference on Structural Mechanics in Reactor Technology (SMiRT15), Seoul, Korea, August 1999, VIII, 13–18.
Dalguer, L. A., K. Irikura, J. D. Riera, and H. C. Chiu (2001). Fault dynamic
rupture simulation of the hypocenter area of the thrust fault of the1999 Chi-Chi (Taiwan) Earthquake, Geophys. Res. Lett. 28, 1327–
1330.
Das, S., and K. Aki (1977). A numerical study of two-dimensional spontaneous rupture propagation, Geophys. J. R. Astr. Soc. 50, 643–668.
Day, S. M. (1982a). Three-dimensional finite difference simulation of fault
dynamics: rectangular fault with fixed rupture velocity, Bull. Seism.
Soc. Am. 72, 705–727.
Day, S. M. (1982b). Three-dimensional simulation of spontaneous rupture:
the effect of nonuniform prestress, Bull. Seism. Soc. Am. 72, 1881–
1902.
Dieterich, J. H. (1979). Modeling of Rock Friction, Part I. Experimental
Results and Constitutive Equation, J. Geophys. Res. 84, 2161–2168.
Doz, G. N., and J. D. Riera (1995). Towards the numerical simulation ofseismic excitation, Transaction, 13th International Conference on
Structural Mechanics in Reactor Technology (SMiRT13), Porto Alegre, Brazil, August 1995, III.
Fukuyama, E., and R. Madariaga (1998). Rupture dynamics of a planar
fault in a 3D elastic medium: rate- and slip-weakening friction, Bull.
Seism. Soc. Am. 88, 1–17.
Guatteri, M., and Spudich, P. (2000). What can strong-motion data tell us
about slip-weakening fault-friction laws? Bull. Seism. Soc. Am. 90,
98–116.
Harris, R. A., and S. M. Day (1999). Dynamic 3D simulations of earthquakes on en echelon faults, Geophys. Res. Lett. 26, 2089–2092.
Ida, Y. (1972). Cohesive force across the tip of a longitudinal-shear crack
and Griffith’s specific surface energy, J. Geophys. Res. 77, 3796–
3805.
Ide, S., and Takeo M. (1997). Determination of constitutive relations offault slip based on seismic wave analysis, J. Geophys. Res. 102,27379–27391.
Inoue, T., and T. Miyatake (1998). 3D simulation of near-field strong
ground motion based on dynamic modeling, Bull. Seism. Soc. Am. 88,
1445–1456.
Iwata, T., H. Sekiguchi, and K. Irikura (2000). Source process of the 1999
Chi-Chi, Taiwan, earthquake and its near-fault strong ground motion,2000 Western Pacific Geophysics Meeting (supplement) EOS Trans.
Am. Geophys. Union 81, WP144, S51F-02 (abstract).
Kostrov, B. V. (1964). Self-similar problem of propagation of shear crack,
J. Appl. Math. Mech. 28, 1077–1087.
Kostrov, B. V. (1966). Unsteady propagation of longitudinal shear cracks,
J. Appl. Math. Mech. 30, 1241–1248.
Ma, K. F., and J. Mori (2000). Rupture process of the 1999 Chi-Chi, Taiwan, earthquake from direct observations and joint inversion of strong
motion, GPS and teleseismic data, 2000 Western Pacific Geophysics
Meeting, Tokyo, Japan, 27–30 June 2000 (supplement) (abstract) EOS
Trans. Am. Geophys. Union 81, WP104, S21A-07.
Madariaga, R., K. Olsen, and R. Archuleta (1998). Modeling dynamic rupture in a 3D earthquake fault model, Bull. Seism. Soc. Am. 88, 1182–
1197.
Marone, C., and Kilgore, B. (1993). Scaling of the critical slip distance forseismic faulting with shear strain in fault zone, Nature 362, 618–621.
Mikumo, T., and T. Miyatake (1978). Dynamical rupture process on a threedimensional fault with non-uniform frictions and near-field seismic
waves, Geophys. J. R. Astr. Soc. 54, 417–438.
Mikumo, T., and T. Miyatake (1993). Dynamic rupture processes on dipping fault, and estimates of stress drop and strength excess from theresults of waveform inversion, Geophys. J. Int. 112, 481–496.
Mora, P., and D. Place (1994). Simulation of the frictional stick-slip instability, PAGEOPH 143, 61–87.
Nayfeh, A. H., and M. S. Hefsy (1978). Continuum modeling of threedimensional truus-like space structures, AIAA J. 16, 779–787.

(A1)

Nielsen, S. B. (1998). Free surface effects on the propagation of dynamic
rupture, Geophys. Res. Lett. 25, 125–128.
Oglesby, D. D., R. J. Archuleta, and S. B. Nielsen (1998). Earthquakes on
dipping faults: the effects of broken symmetry, Science 280, 1055–
1059.
Oglesby, D. D., R. J. Archuleta, and S. B. Nielsen (2000). The threedimensional dynamics of dipping faults, Bull. Seism. Soc. Am. 90,
616–628.
Ohnaka, M. (1992). Earthquake source nucleation: a physical model for
short term precursors, Tectonophysics 211, 149–178.
Ohnaka, M., Y. Kuwahara, and K. Yamamoto (1987). Constitutive relations
between dynamic physical parameters near a tip of the propagating
slip zone during stick-slip shear failure, Tectonophysics 144, 109–
125.
Olsen, K., R. Madariaga, and R. Archuleta (1997). Three-dimensional dynamic simulation of the 1992 Landers earthquake, Science 278, 834–
838.
Riera, J. D. (1982) Basic concepts and load characteristics in impact problems, in Concrete Structures under Impact and Impulsive Loading,Introductory report, RILEM, CEB, IABSE, IASS-Interassociation
Symposium, Berlin, June 1982, 7–29.
Riera, J. D., and M. Rocha (1991). A note on the velocity of crack propagation in tensile fracture, Revista Brasileira de Ciencias Mecanicas,
RBCM, XII, no. 3, 217–240.
Rimal, M. H. (1992). Dynamic fracture analyses by the extended distinct
element Method, Doctoral Dissertation, Department of Civil Engineering, University of Tokyo, Japan.
Ruina, A. (1983). Slip instability and state variable friction laws, J. Geophys. Res. 88, 10359–10370.
Scholz, C. H. (1990). The Mechanics of Earthquakes and Faulting, Cambridge University Press, New York.
Shi, B., A. Anooshehpoor, J. N. Brune, and Y. Zeng (1998). Dynamics
thrust faulting: 2D lattice model, Bull. Seism. Soc. Am. 88, 1484–
1494.
Shin, T. C., K. W. Kuo, W. H. K. Lee, T. L. Teng, and Y. B. Tsai (2000).
A preliminary report on the 1999 Chi-Chi (Taiwan) Earthquake,
Seism. Res. Lett. 71, 24–30.
Tsai, Y. B., and M. W. Huang (2000). Strong ground motion characteristicsof the Chi-Chi, Taiwan earthquake of September 21, 1999, Earthquake Eng. Eng. Seismology 2, 1–21.
Virieux, J., and R. Madariaga (1982). Dynamic faulting studied by a finite
difference method, Bull. Seism. Soc. Am. 72, 345–369.
Yagi, Y., and M. Kikuchi (2000). Source rupture process of the Chi-Chi,
Taiwan, earthquake of 1999, obtained by seismic wave and GPS data,
2000 Western Pacific Geophysics Meeting, Tokyo, Japan, 27–30 June
2000 (supplement) (abstract), LOS Trans. Am. Geophys. Union 81,
WP104, S21A-05.

Appendix
Determination of the Equivalent Stiffness of the Normal
and Diagonal Bars of the DEM
The determination of the equivalent axial stiffness ofthe elements in a cubic lattice array, as shown in Figure 4(equations 1 and 2), in terms of the elastic properties of anequivalent isotropic continuum is reviewed here. This equivalence was shown by Nayfeh and Hefzy (1978) and first
employed in dynamic problems by Riera (1982).
The stress-strain for a general elastic body may be written in the indicial notation form:ri ⳱ Cij ej , i,j ⳱ 1 . . . 6,

1126

L. A. Dalguer, K. Irikura, J. D. Riera, and H. C. Chiu

where ri and ej are the independent six components of thestress and strain tensors, respectively, and Cij are elastic constants, which are 21 independent constants on account ofsymmetry considerations.
For an isotropic material, the constants Cij could be afunction of only two independent constants, so Cij can beordered as:

Cij ⳱

冤

C11 C12 C12 0
00
C12 C11 C12 0
00
C12 C12 C11 000
0
0
0 C44 000
0
0
0 C44 0
0
0
0
0
0 C44

冥

,

(A2)

in which C11, C12, and C44 are functions of Young’s modulus
E and the Poisson’s ratio v.
Since the DEM implies a lattice-type structure consisting of one-dimensional axial elements, the contribution ofeach member to the overall stiffness must be duly accounted
for the sum of the average contribution of each element willlead to the final stiffness matrix. It is assumed that theelements are perfectly straight and present constant crosssectional area.
The elastic constants Cij can be transformed from one
orthogonal cartesian coordinates system xi to another x̄i (i ⳱
1,2,3) through an expression:
Q̄ij ⳱ f(Qi,j, ␣kl)冦(k, l ⳱ 1 . . . 3a),(i, j ⳱ 1 . . . 6)

n
Q11
⳱

d
Q11
⳱

2EAn
,
L2

冪3EAd
L2

(A5)

,

in which L is the side length of the basic cube, and EAn andEAd are Young’s modulus E times the cross-sectional area
A for the normal and diagonal elements, respectively.
From equations (A4) and (A5), Q̄ij of equation (A3),that is, the stiffness constants of a continuum equivalent tothe cubic lattice model may be obtained. As shown in Figure
4a, seven elements converge to each node of the cubic model(three normal and four diagonals bars), hence Q̄ij can be
expressed as follows:3

Q̄ij ⳱

f(Qn11, ␣nIkl)
兺
I⳱1
4

Ⳮ

兺 f(Qd11, ␣dJkl) (k, l ⳱ 1 . . . 3),

(A6)

J⳱1

in which ␣nIkl and ␣dJkl are the direction cosines of the systems
x̄, xnI and x̄, xdJ , respectively.
Combining equations (A4), (A5), and (A6), the elastic
constants of equation (A2) can be obtained:

冢

冣

4
d ,
9

n
C11 ⳱ Q11
1 Ⳮ

(A3)

where ␣n denotes the direction cosines of the transformation.
Qij and Q̄ij are the elastic constant of systems x and x̄, respectively. The derivation of equation (A3) is presented indetails in Nayfeh and Hefzy (1978).
Equation (A3) is used in order to get the coefficients ofequation (A2) for the corresponding cubic model of Figure
4a. Since all elements have the single unidirectional property
E, each set of parallel bars defines a continuum with a single
effective unidirectional property, which we shall refer to as
Q11. In the context of effective modulus theories, Q11 will
be an area-averaged modulus. Thus, the value of Q11 will
depend not only upon the specific model under consideration
but also upon the spacing of the bars.
In particular, the cubic array shown in Figure 4a has
two different properties of Q11, one corresponds to the elements normal to the face of the cube (Qn11), and the othercorresponds to the diagonal elements (Qd11).
The unidirectional effective properties, Qn11 and Qd11, canbe determined by referring to the projected area normal toeach element in order to obtain the effective area. For thenormal elements, we can obtain an effective area equal toL2/2, while for the diagonal elements it results in L2/冪3.
Thus, the unidirectional property for each element is given
by

(A4)

n
C12 ⳱ Q11

冢94 d冣,n
C44 ⳱ Q11

冢49 d冣,

(A7)

where
d⳱

冪3Ad
Qd11
.
n ⳱
Q11
2An

(A8)

Therefore, the constants Cij of equation (A2) is:

冤

4d1 Ⳮ 㛭9

Cij ⳱

2EAn
L2

4d
㛭
9
4d
㛭
9

4d
㛭

4d
㛭

9
4d
㛭

1Ⳮ 9
4d
㛭
9

9
4d
㛭
9
4d
㛭

1Ⳮ 9

0

0

0

0

0

0

0

0

0

0

0
0

0

0

0

4d
㛭

0

0

0

0

4d
㛭

0

0

0

0

0

9

9

4d
㛭
9

冥

(A9)

The engineering elastic constants, Young’s modulus E,
Poisson’s ratio v, and shear modulus l may be expressed interms of the constants Cij using the general stress-strain re-

The Importance of the Dynamic Source Effects on Strong Ground Motion during the 1999 Chi-Chi, Taiwan, Earthquake

lations for anisotropic materials (Nayfeh and Hefzy, 1978),as follows:
E ⳱ C11 ⳮ

m⳱

2
2C12
,
C11 Ⳮ C12

(A10a)

C12,
C11 Ⳮ C12

(A10b)1
C44 .
2(A10c)l⳱EAn
L(9 Ⳮ 8d)⳱ELn2(9 Ⳮ 12d)

(A13a)EAd
2 dL(9 Ⳮ 8 d)⳱ELd3(9 Ⳮ 12 d)(A13b)The nodal masses, mi, are calculated in terms of thevolumes of influence of each node i, for the internal nodes
it results:mi ⳱

Substituting in equation (A9), we obtain:
2EAn (9 Ⳮ 12d),L2(9 Ⳮ 8d)(A11a)m⳱

4d,
9 Ⳮ 8d(A11b)l⳱

EAn 4d.
9L2(A11c)

E⳱

To get the equivalent values of EAn and EAd as functionsof the elastic properties of the material (which could be, forinstance, E and v), we use equation (A11) and the relation
given in equation (A8), as follows:
9m,(4 ⳮ 8m)(A12a)

L2(9 Ⳮ 8d)E,2(9 Ⳮ 12d)

(A12b)

d⳱

EAn ⳱

EAd ⳱

dL2(9 Ⳮ 8d)

冪3(9 Ⳮ 12d)

E.

(A12c)

Finally, the equivalent stiffness of the normal and diagonal elements is obtained dividing equations (A12a,b) by
the length of the respective element, being Ln ⳱ L (for normal elements) and Ld ⳱ 冪3/2 L (for diagonal elements):

1127

qL3
2

(A14)

while for surface nodes, linear corner, and point corner
nodes, the mass given by equation (A14) must be divided
by 2, 4, and 8, respectively.
Disaster Prevention Research InstituteKyoto University
Gokasho, Uji, Kyoto 611-0011, Japan
dalguer@egmdpri01.dpri.kyoto-u.ac.ip(L.A.D.)
Curso de Pós Graduação de Engenharia Civil
Universidade Federal de Rio Grande do Sul
Av. Osvaldo Aranha 99, 3⬚ andar, CEP 90035-190
Porto Alegre, R.S. Brazil
dalguer@vortex.ufrgs.br(L.A.D.)
Disaster Prevention Research InstituteKyoto University
Gokasho, Uji, Kyoto 611-0011, Japan
irikura@egmdpri01.dpri.kyoto-u.ac.ip(K.I.)
Curso de Pós Graduação de Engenharia Civil
Universidade Federal de Rio Grande do Sul
Av. Osvaldo Aranha 99, 3⬚ andar, CEP 90035-190
Porto Alegre, R.S. Brazil
riera@genesis.cpgec.ufrgs.br(J.D.R.)
Institute of Earth Science
Academia Sinica, Nankang, Taipei, Taiwan
chiu@earth.sinica.edu.tw(H.C.C.)
Manuscript received 31 July 2000.