of 32
1765
Bulletin of the Seismological Society of America, 91, 6, pp. 1765–1796, December 2001
Dynamic Earthquake Ruptures in the Presence of Lithostatic Normal
Stresses: Implications for Friction Models and Heat Production
by Brad T. Aagaard, Thomas H. Heaton, and John F. Hall
Abstract
We simulate dynamic ruptures on a strike-slip fault in homogeneous
and layered half-spaces and on a thrust fault in a layered half-space. With traditional
friction models, sliding friction exceeds 50% of the fault normal compressive stress,
and unless the pore pressures approach the lithostatic stress, the rupture character-
istics depend strongly on the depth, and sliding generates large amounts of heat.
Under application of reasonable stress distributions with depth, variation of the ef-
fective coefficient of friction with the square root of the shear modulus and the inverse
of the depth creates distributions of stress drop and fracture energy that produce
realistic rupture behavior. This
ad hoc
friction model results in (1) low-sliding fric-
tion at all depths and (2) fracture energy that is relatively independent of depth.
Additionally, friction models with rate-weakening behavior (which form pulselike
ruptures) appear to generate heterogeneity in the distributions of final slip and shear
stress more effectively than those without such behavior (which form cracklike rup-
tures). For surface rupture on a thrust fault, the simple slip-weakening friction model,
which lacks rate-weakening behavior, accentuates the dynamic interactions between
the seismic waves and the rupture and leads to excessively large ground motions on
the hanging wall. Waveforms below the center of the fault (which are associated
with waves radiated to teleseismic distances) indicate that source inversions of thrust
events may slightly underestimate the slip at shallow depths.
Introduction
A better understanding of the rupture process provides
an important avenue for improving the models of near-
source ground motions and gaining insight into the physics
of earthquakes. Including the rupture dynamics in simula-
tions of earthquakes generally involves modeling the fric-
tional sliding on the fault surface, with two distinct efforts
having emerged in recent years. Those researchers who
model the evolution of stress on the fault almost exclusively
use state- and rate-dependent friction models. Review arti-
cles by Marone (1998) and Scholz (1998) summarize the
development of the friction models and some of the features
of their behavior. These models are based on laboratory ex-
periments of sliding at slip rates between 10

7
and 1 mm/
sec and can be derived from analytical models of creep be-
havior (Persson, 1997). Consequently, researchers apply
these models to studies of the nucleation of earthquakes and
creep behavior on faults (e.g., Stuart and Mavko, 1979; Rice
and Ben-Zion, 1996; Tullis, 1996). As the coefficient of fric-
tion in these state- and rate-dependent friction models is typ-
ically around 0.6, their application to dynamic ruptures pre-
dicts large temperature changes in the zone surrounding the
fault, unless the dynamic compressive stresses on the fault
are much less than the lithostatic pressures (Richards, 1976;
Kanamori and Heaton, 2000).
The other effort focuses on modeling the rupture be-
havior during earthquakes. The uncertainty in the behavior
of how faults rupture has led researchers to create simple,
ad hoc
models that produce reasonable behavior. These
models generally include either slip-weakening behavior
(the shear strength decreases as slip occurs) or a combination
of slip weakening and rate weakening (initially, the shear
strength drops with slip in response to slip weakening and
then returns near its original level as the slip rate decreases).
For nearly 30 yr, the slip-weakening friction model has been
used to study the frictional sliding associated with earth-
quakes. Ida (1972) was one of the first to associate slip-
weakening behavior with the propagation of shear cracks.
Andrews (1976a) and Burridge
et al.
(1979) used slip-
weakening friction models to study the propagation of mode-
II shear cracks.
At about the same time, the finite-difference and finite-
element methods were applied to the study of three-dimen-
sional dynamic ruptures (Madariaga, 1976; Archuleta and
Day, 1980; Day, 1982a). This marked a dramatic improve-
1766
B. T. Aagaard, T. H. Heaton, and J. F. Hall
ment in the applicability of the methods used to model dy-
namic ruptures because they can be used for three-dimen-
sional simulations with heterogeneous material properties;
however, the computing power at the time severely limited
the size and scope of the calculations.
More recently, with the advances in computing, many
more researchers have employed boundary-integral, finite-
difference, or finite-element formulations to model dynamic
ruptures (Mikumo, 1992; Mikumo and Miyatake, 1993; Ma-
dariaga and Cochard, 1996; Harris and Day, 1997; Olsen
et
al.,
1997; Fukuyama and Madariaga, 1998; Oglesby
et al.,
1998; Harris and Day, 1999; Magistrale and Day, 1999;
Nielsen and Olsen, 2000; Oglesby
et al.,
2000a,b). Simula-
tions of the 1992 Landers earthquake by Olsen
et al.
(1997)
and the 1994 Northridge earthquake by Nielsen and Olsen
(2000) demonstrated the ability of the finite-difference
method and slip-weakening friction models and slip- and rate-
weakening friction models, respectively, to produce reason-
able rupture behavior. The simulations generated a confined
rupture pulse consistent with the kinematic source inversions
and reproduced the main long-period features of the wave-
forms. Using similar methods, Harris and Day (1999) and
Magistrale and Day (1999) explained the propagation across
step-overs between parallel strands on strike-slip and thrust
faults, respectively. Additionally, Oglesby
et al.
studied the
difference between ruptures on normal faults and thrust
faults using three-dimensional (Oglesby
et al.,
2000a,b)
finite-element simulations in a homogeneous medium.
In contrast to the simulations using state- and rate-
dependent friction models, only a few of the simulations
performed to date with slip-weakening or slip- and rate-
weakening friction models assume effective normal stresses
that correspond to the overburden pressure. Although Mik-
umo (1992) used normal stresses equal to the overburden
pressure, the distributions of final slip from the dynamic rup-
tures exhibit a clear depth dependence, which does not match
the distributions found in kinematic source models (Heaton,
1990; Somerville
et al.,
1999). Consequently, we examine
what constraints on the friction model may be required to
produce realistic ruptures when we apply reasonable normal
stress distributions with depth. Using this implementation of
the earthquake source and following the constraints imposed
on the friction model, we determine the sensitivity of the
rupture behavior and the ground motions to systematic var-
iations of the initial shear stresses and the friction model.
Additionally, we compare the ground motions from the dy-
namic rupture simulations with the corresponding cases of
prescribed ruptures. Aagaard (1999) includes an expanded
discussion of dynamic ruptures in both homogeneous and
layered half-spaces.
Methodology
We use the same general solution techniques, which are
described by Aagaard (1999) and Aagaard
et al.
(2001), to
simulate the earthquakes with dynamic ruptures that we use
for prescribed ruptures. Applying the finite-element method
with linear tetrahedral elements to the three-dimensional dy-
namic elasticity equation produces the matrix differential
equation,
[ ]{ ( )} [ ]{ ( )} [ ]{ ( )} { ( )},
Mu
Cu
Ku
F

tt tt
++=
(1)
where [
M
], [
C
], and [
K
] denote the mass, damping, and
stiffness matrices, respectively, {
u
(
t
)} denotes the displace-
ment vector at time
t
, and {
F
(
t
)} denotes the force vector at
time
t.
We model the slip on the fault using sliding degrees
of freedom (see Aagaard, 1999; Aagaard
et al.,
2001 for
details) which create dislocations on the fault surface. On
the fault surface we transform the usual three translational
degrees of freedom on each side of the fault to six fault
degrees of freedom consisting of three relative degrees of
freedom (two tangential to the fault plane and one normal
to the fault plane) and three average degrees of freedom (two
tangential to the fault plane and one normal to the fault
plane). This transformation provides explicit control of the
relative motion between the two sides of the fault. In pre-
scribed ruptures the specified slip time histories dictate the
relative tangential displacements of the fault degrees of free-
dom, whereas in dynamic ruptures the friction model con-
strains the forces acting on the relative tangential fault de-
grees of freedom. Incorporating the fault surface into the
geometry of the finite-element model allows arbitrary ori-
entation of the fault plane.
The seismic waves generated by the rupturing fault
create dynamic stresses in the surrounding volume as well
as changes in the static stresses. We assume that the static
stresses on the boundary of the domain remain constant dur-
ing the earthquake (see the Appendix for a discussion on
how this affects the energy balance). We do not need to
know the initial stresses throughout the domain to model the
seismic-wave propagation. However, in order to simulate the
dynamic rupture of the fault, we must know the initial
stresses acting on the fault surface. These stresses may be
found in a number of ways, including solution of a static
problem, solution of a viscoelastic problem, extrapolated
from data, or assumed from intuition. Regardless of their
source, we resolve the stresses into shear and normal trac-
tions acting on the fault surface. Thus, off the fault surface
we consider only the dynamic stresses and the change in the
static stresses, whereas on the fault surface we also consider
the initial (static) stresses.
We treat the friction on the fault as an external force
and replace the force vector in the governing equation (equa-
tion 1) with the difference between the vector of tectonic
forces, {
F
t
}, and the friction force vector, {
F
f
}; this yields
[ ]{ ( )} [ ]{ ( )} [ ]{ ( )} { ( )} { ( ( ), ( ))},
Mu
Cu
Ku
F
F


t
t
t
t
Dt Dt
++=−
tf
(2)
where
D
denotes the slip on the fault. As outlined previously,
we only apply the tectonic forces to the fault degrees of
freedom. The appearance of the difference between the tec-
Dynamic Earthquake Ruptures in the Presence of Lithostatic Normal Stresses
1767
tonic force vector and the friction force vector in the equa-
tion of motion implies that we may create the same sliding
behavior from an infinite combination of tectonic and fric-
tion forces by keeping the difference between them the same.
The vector {
F
f
} acts on the relative tangential fault de-
grees of freedom, whereas the vector {
F
t
} acts on both the
relative normal and relative tangential fault degrees of free-
dom. Following the same procedure that we use for pre-
scribed ruptures, we integrate the differential equation using
the central-difference scheme. When the coefficient of fric-
tion depends on the slip rate, computing the friction at time
t
requires knowing the slip rate at time
t
, which we do not
know. To remedy this difficulty, we assume that the time
step is small enough so that the slip rate does not change
significantly in a single time step. Thus, we use the slip rate
at time
t

D
t
, instead of the slip rate at time
t
, to compute
the friction force at time
t.
We must transform the tectonic tractions applied on the
fault surface into forces acting on the fault degrees of free-
dom. At each node on the fault, we interpolate from a given
distribution of initial tractions and convert the tractions to
forces using the node’s tributary area on the fault plane. We
assume that the fault is in equilibrium and apply the forces
to the relative degrees of freedom. The friction force does
not require any transformation; the product of the coefficient
of friction and the force associated with the relative normal
degree of freedom gives the maximum magnitude of the fric-
tion force vector acting on the relative tangential degrees of
freedom. The dynamic deformation in the domain may cause
variations in the normal forces acting on the fault, but we
do not allow normal separation of the two sides of the fault.
Furthermore, except at the ground surface, the confining
pressures keep the normal stresses well within the compres-
sive regime. We compute the dynamic normal forces at the
fault degrees of freedom as part of the time-stepping pro-
cedure.
Initial Tractions on Fault
We consider gravity and plate tectonics as sources of
normal stresses acting on the fault surface. In a self-gravi-
tating, spherical Earth with only radial variations in material
properties, the weight of the material generates lithostatic
stresses (total stress caused by gravity) with no shear stresses
and equal axial stresses (Mohr’s circle degenerates into a
point). In addition to shear stresses, plate tectonics also cre-
ates normal stresses on the fault surface, especially in the
case of inclined faults.
The presence of water in the interstices of the granular
medium can generate pore pressures that decrease the effec-
tive normal stresses. If little or no water sits in the interstices,
then the effect of the pore pressures is negligible, so the
effective normal stresses equal the normal stresses. In a dry,
homogeneous half-space, the effective normal stresses (lith-
ostatic stresses) increase linearly with depth (
p

q
gz
). If
water saturates the interstices, then the pore pressures equal
the hydrostatic pressures, and the effective normal stresses
are the difference between the normal stresses and the hy-
drostatic pressures. In a saturated, permeable, homogeneous
half-space, the effective normal stresses again increase lin-
early with depth, but at a slower rate because of the presence
of hydrostatic pore pressures [
p

(
q

q
w
)
gz
]. Finally, if
the rock is saturated but impermeable, the pore pressures can
equal the lithostatic pressures, and the effective normal
stresses can become very small; in this case the material
essentially floats. The existence of topography and density
variations implies large shear stresses (

10 MPa) at depths
that require large normal stresses to prevent failure. Conse-
quently, except in localized areas, we expect the effective
normal stresses to be similarly large. Researchers often use
effective normal stresses that are independent of depth for
simplicity (Olsen
et al.,
1997; Ben-Zion and Andrews, 1998;
Madariaga
et al.,
1998), without acknowledging that assum-
ing uniform effective normal stresses with depth implies that
the pore pressures increase more rapidly than the hydrostatic
pressures such that the differences between the lithostatic
pressures and the pore pressures are uniform with depth.
Shear tractions on the fault generate the forces that
cause slip on the fault surface. We apply the shear tractions
to the relative tangential degrees of freedom in the direction
of the desired slip and use an asperity (usually circular in
shape) with a shear stress greater than the failure stress to
start the rupture. Many factors, such as the discretization
size, the failure stress, and the dynamic stress drop, influence
the size of the asperity necessary to initiate a propagating
rupture (Andrews, 1976b; Day, 1982b; Madariaga
et al.,
1998).
Overview of Rupture Dynamics
We examine the anatomy of the shear stress on the fault
near the rupture front, as shown in Figure 1, to find the
relationship between its features and the dynamics of the
rupture. From an elasticity theory treatment of dynamic frac-
ture mechanics (see Freund, 1990), the shear stresses in-
crease and form a singularity just ahead of the leading edge
of the rupture. After slip begins, the shear stresses decrease,
and then, depending on the friction model, may or may not
recover as the slip rate decreases. However, in our finite-
element models (as in all discrete models) the shear stresses
remain finite with a stress concentration at the rupture front
rather than a singularity.
The friction model controls the decrease in friction
stress as slip progresses, and therefore the dynamic stress
drop (the difference between the initial shear stress and the
shear stress during sliding) as well. The rate at which the
dynamic stress drop increases behind the rupture front gov-
erns the slip rate, with faster decreases in shear stress leading
to greater slip rates. Additionally, a larger dynamic stress
drop results in a larger stress concentration at the leading
edge of the rupture. The increase in shear stress associated
with the stress concentration dictates when slip occurs at
each point and, as a result, the rupture speed. Thus, the rup-
1768
B. T. Aagaard, T. H. Heaton, and J. F. Hall
from failure
"Distance"
Stress
drop
Dynamic
stress
drop
Initial
stress
Slip ends
Slipping region
Stress
Position
Direction of propagation
Stress
Slip begins
Time
Figure
1.
Diagrams of the concentration of shear stress near the rupture front at a
specific time as a function of space (left diagram) and at a fixed location as a function
of time (right diagram).
D
σ
f
σ
b
fa
σ
m
σ
a
fa
D
b
o
D
a
o
Figure
2.
Friction stress (
r
f
) as a function of slip
distance (
D
) with an illustration of two sets of param-
eters for the slip-weakening friction model (denoted
by the superscripts a and b) that have the same frac-
ture energy (hatched areas), but different failures
stresses (
and
) and slip distances (
and
).
ab
ab
rr
DD
fail
fail
0
0
ture speed and slip rate are not independent but are related
through the dynamic stress drop.
We may also consider the dynamics of the rupture using
energy. As the rupture propagates, the rupture front con-
sumes energy through sliding. We associate two forms of
energy with the sliding. We call the energy dissipated when
the friction decreases during sliding the fracture energy (il-
lustrated in Fig. 2) because it corresponds to the fracture
energy in crack models (Rice, 1983). We associate the en-
ergy dissipated through sliding at a relatively constant fric-
tion stress with the generation of heat. The sliding also re-
leases the energy radiated in the seismic waves. As we
increase the fracture energy for a given maximum dynamic
stress drop, the rupture consumes more energy leaving less
available for radiation. In such cases the slip rates and rup-
ture speed decrease (Rice, 1983; Fukuyama and Madariaga,
1998). Likewise, with a decrease in the fracture energy, more
energy is available for sliding, and the slip rates and rupture
speed increase. If the fracturing dissipates more energy than
the energy released, then the rupture slows and eventually
stops. The rupture will also stop propagating if the leading
edge of the rupture slows down and it is caught from behind
by the trailing edge of the rupture.
Conceptually, we wish to separate the energy dissipated
on the fault surface into fracture energy and heat generated
by the frictional sliding. However, most conventional fric-
tion models define a slip-weakening distance
D
0
, where the
shear stress decreases from the failure stress to the frictional
sliding stress over the slip distance
D
0
. If the shear stress
decreases linearly over the slip-weakening distance, the frac-
ture energy per unit area equals one-half of the strength ex-
cess times the slip-weakening distance (see Fig. 2). In dis-
crete models the slip-weakening distance cannot become
arbitrarily small because of the finite discretization size.
Figure 2 shows that the failure stress does not uniquely
determine the fracture energy: adjusting the slope of the fric-
tion model changes the fracture energy. For example, at a
lower failure stress, we can maintain the same fracture en-
ergy by reducing the rate at which the friction stress de-
creases with slip. Guatteri and Spudich (2000) demonstrated
this method for an event resembling the
M
6.5 1979 Imperial
Valley, California, earthquake.
This technique plays a critical role in manipulating the
dynamics of the rupture in simulations with discretized do-
mains. Accurately capturing the stress concentration in shear
stress near the leading edge of the rupture requires much
finer discretization than that necessary to model the wave
propagation (Day, 1982b; Madariaga
et al.,
1998) because
the failure stress develops only over a very localized region.
We wish to capture the general features of such failure with-
out modeling such localized behavior. In a discrete model,
such as a finite-element model, the failure stress becomes a
parameter dependent on discretization size, but the fracture
energy should continue to control the behavior of the rup-
Dynamic Earthquake Ruptures in the Presence of Lithostatic Normal Stresses
1769
0
1
D/D
o
Coefficient of Friction
μ
min
μ
max
fracture
energy
Figure
3.
Slip-weakening friction model. The co-
efficient of friction decreases over the slip distance
D
0
. The shaded region is associated with the fracture
energy.
Table 1
Description of the Variables Involved in the Friction Models
Variable
Dimensions
Description
l
f
Dimensionless
Coefficient of friction
l
max
Dimensionless
Maximum coefficient of friction
l
min
Dimensionless
Minimum coefficient of friction
D
Length
Slip distance
V
,
̇
D
Length/time
Slip rate
D
0
Length
Slip distance constant
V
0
Length/time
Slip rate constant
ture. We can manipulate the friction model to maintain the
same fracture energy for different levels of the failure stress
by changing the slip distance as demonstrated in Figure 2.
Clearly, this technique breaks down when the slip-weaken-
ing distance,
D
0
, exceeds the actual amount of slip. Addi-
tionally, for a given fracture energy, as the slip-weakening
distance and discretization size increase, the rupture jumps
more easily to super-shear rupture speeds. Despite these dif-
ficulties, we can use larger finite elements than those re-
quired to accurately capture the stress concentration and al-
low the wave propagation to control the local element size
without significantly altering the behavior of the rupture. In
other words, the failure shear stress in our finite-element
models, which determines when a point on the fault begins
to slip, is actually some measure of the shear stress at failure
averaged over the discretization size. Consequently, it does
not correspond to the yield stress in a continuum.
Friction Models
We focus on two models of sliding friction, both of
which compute the upper bound on the friction force from
the product of the normal force and the coefficient of fric-
tion. We do not implement the state- and rate-dependent
friction models advocated by several researchers, including
Dieterich (1992) and Scholz (1998), because they are ob-
served at slip rates of less than 1 mm/sec and are associated
with viscoelastic creep behavior (Persson, 1997); slip during
earthquakes occurs at rates of the order of meters per second
(Heaton, 1990; Somerville
et al.,
1999). Furthermore, state-
and rate-dependent friction models imply high-sliding fric-
tion in the presence of large fault normal compressive
stresses. Several plausible mechanisms have been suggested
to explain why friction varies in both space and time during
earthquakes for dynamic reasons, for example, wrinklelike
slip pulses associated with a contrast in material properties
(Harris and Day, 1997; Ben-Zion and Andrews, 1998; An-
ooshehpoor and Brune, 1999), acoustic fluidization (Melosh,
1996), normal vibrations (Brune
et al.,
1993; Tworzydlo and
Hamzeh, 1997), and elastohydrodynamic lubrication (Brod-
sky and Kanamori, 2001). Instead of choosing to model any
particular mechanism, we approximate the general features
of the behaviors with simple friction models because several
of these mechanisms may be combining to change the fric-
tion stress during sliding. Thus, we choose to use simple,
ad
hoc
friction models with characteristics such as slip weak-
ening (decrease in friction with the progression of slip) and
rate weakening (increase in friction as the slip rate ap-
proaches zero) that produce realistic rupture behavior and
capture the general features of the more complicated models.
Table 1 provides descriptions of the parameters used in the
functional forms of the coefficient of friction.
The two end member cases of rupture models include
cracklike behavior, where the healing phases emanate from
the boundaries of the fault, and pulselike behavior, where
healing phases occur spontaneously and trail behind the
leading edge of the rupture (Heaton, 1990). In slip-weak-
ening friction models that model true cracklike behavior, the
sliding friction drops to some level as slip occurs and re-
mains there; no restrengthening occurs even when sliding
stops. With no shear restrengthening, the slip tends to over-
shoot its final value, and slip occurs whenever seismic waves
with a nonzero fault tangential component attempt to prop-
agate across the fault. The slip can by reduced to generally
only one episode during a rupture by including shear res-
trengthening, wherein the sliding friction returns to its initial
value upon termination of sliding. As we shall see, pulselike
behavior, instead of cracklike behavior, develops when we
allow a more gradual increase in friction as the slip rate tends
toward zero (rate-weakening behavior).
In the slip-weakening friction model the coefficient of
friction decreases linearly from a maximum value to a min-
imum value over a slip distance of
D
0
:
μ
μ
μμμ
μ
f
max
max
max
min
min
=
=
−−≤
>

Dt
Dt
D
Dt
D
Dt
D
()
()
()
()
()
.
0
0
0
0
(3)
This defines the latent heat (fracture energy) generated by
fracture for this friction model. We refer to this model as
slip-weakening friction because the material exhibits a
weakening in shear strength as slip occurs. Figure 3 illus-
trates how the coefficient of friction decreases from
l
max
to
l
min
over a slip distance of
D
0
. When the slip rate returns to
1770
B. T. Aagaard, T. H. Heaton, and J. F. Hall
0
1
2
3
4
0
1
2
3
4
D/D
o
V/V
o
μ
max
μ
post
μ
min
Coefficient of Friction
Figure
4.
Slip- and rate-weakening friction model.
The thick line indicates a typical trajectory of the co-
efficient of friction. The coefficient of friction de-
creases over the slip distance
D
0
and then increases
when the slip rate drops below
V
0
. The shaded region
is associated with the fracture energy.
zero, we allow shear restrengthening, so the coefficient of
friction returns to
l
max
. This results in a smaller static stress
drop compared with the dynamic stress drop. Without shear
restrengthening the static stress drop may exceed the dy-
namic stress drop (Madariaga, 1976).
Following Madariaga and Cochard (1996) and Mada-
riaga
et al.
(1998), we create a second friction model that
depends on slip distance and slip rate by taking the greater
of the two coefficients of friction determined from the slip-
weakening friction model and a rate-weakening friction
model. As a result, there is no simple expression for the
coefficient of friction as a function of slip distance and slip
rate. The rate-weakening friction model corresponds to re-
placing the slip distance in the slip-weakening friction model
with the slip rate and the slip distance
D
0
with the slip rate
V
0
. We also replace
l
max
in the rate-weakening friction
model with
l
post
to allow different shear strengths before and
after slip. We refer to this model as slip and rate weakening.
Figure 4 illustrates the variation of the coefficient of friction
with both slip distance and slip rate, and a typical path during
sliding.
Dynamic Energy Balance
The energy balance provides an additional tool for char-
acterizing an earthquake, and the change in thermal energy
allows estimation of the degree of melting on the fault. We
derive the dynamic energy balance from the conservation of
energy for the entire Earth, assuming no heat is lost on the
timescale of the earthquake. We neglect all external forces,
such as the gravitational forces from the sun and the other
planets, and therefore have no change in the internal energy
of the Earth. As given by
EQW
R
++=
∆∆
0,
(4)
the internal energy of the Earth consists of the radiated en-
ergy (
E
R
), the change in thermal energy (
D
Q
), and the
change in the potential energy (
D
W
). We ignore the rota-
tional energy of the Earth so that the change in potential
energy equals the sum of the change in the strain energy and
the change in the gravitational potential energy.
When we think about energy and earthquakes, we often
only consider the radiated energy because we associate it
with the ground motions and can estimate it from ground-
motion records. Similarly, in our numerical models the ra-
diated energy is readily available from the earthquake simu-
lation by finding the energy dissipated through the damping
matrix.
The primary contribution to the change in the thermal
energy comes from the generation of heat by the frictional
sliding on the fault. Additionally, the fracturing of materials
in the fault zone creates latent heat. The radiated energy
eventually dissipates into heat, but we consider it separately
as discussed earlier. We include both the fracture behavior
and the sliding behavior in the friction model. Consequently,
the energy dissipated through frictional sliding includes both
the latent heat associated with the fracture energy and the
heat generated by sliding. In order to find the energy dissi-
pated during frictional sliding on the fault (
D
Q
(
t
)), we in-
tegrate the increment of heat produced by an increment of
slip over the fault surface:
Qt
tDt S t
S
t
()
() ()
,
=
σ
f
dd

(5)
where
r
f
(
t
) and

Dt
()
are the friction stress and slip rate at
time
t.
The heat generated during sliding will increase the tem-
perature in the region surrounding the fault. Because equa-
tion (5) includes the fracture energy (which does not induce
temperature changes), we use the product of the minimum
friction stress during sliding and the final slip to compute
the heat per unit area (
D
Q
temp
) generated by the sliding at
each point on the fault. We find the change in temperature
at a point on the fault using
T
Q
Cd
=
temp
v
ρ
,
(6)
where
C
v
denotes the heat capacity per unit mass,
q
denotes
the mass density, and
d
is the maximum distance perpendic-
ular to the fault to which the heat penetrates.
The distinction between fracture energy and frictional
heat is important because only the frictional heat increases
the temperature in the sliding zone of the fault. Estimates of
the degree of melting in the sliding zone constrain the
amount of heat that can be generated over the timescale of
the dynamic rupture. On the other hand, fracture energy can
be dissipated by anelastic deformation in the vicinity of the
leading edge of the propagating rupture (crack tip). Esti-
mates of fracture energy in large earthquakes generally fall
in the range of 10
6
J/m
2
, which is 2 to 4 orders of magnitude
larger than the fracture energies observed in the laboratory
Dynamic Earthquake Ruptures in the Presence of Lithostatic Normal Stresses
1771
(Beroza and Spudich, 1988; Lawn, 1993). Such large frac-
ture energies suggest that anelastic yielding occurs over a
much broader region than where the predominant slip takes
place.
We define the change in potential energy as the energy
released by the slip on the fault, assuming that the slip occurs
quasi-statically and that the domain behaves according to
linear elasticity. Because both the radiated energy and the
change in heat energy must be positive, conservation of en-
ergy dictates that the change in potential energy must be
negative. This drop in the potential energy allows earth-
quakes to release energy as propagating waves and to gen-
erate heat through frictional sliding.
We follow a procedure similar to that of Dahlen (1977)
and Savage and Walsh (1978) to find the change in potential
energy caused by an earthquake. Starting with the change in
energy for an increment of slip and assuming that the me-
dium behaves linearly elastically gives us
WDS
S
=−+
1
2
01
(),
σσ
d
(7)
where
D
is the slip, and
r
0
and
r
1
are the shear stresses
before and after the earthquake, respectively.
From the point of view of understanding the physics of
the rupture, we would like to decompose the change in po-
tential energy into change in strain energy and change in
gravitational potential energy. As shown by Dahlen (1977)
and Savage and Walsh (1978), we cannot determine these
changes in energy when we truncate the domain because all
points in the Earth contribute terms of the same order in the
computations; the domain must encompass the entire Earth
in order to compute the change in strain energy or the change
in gravitational potential energy. The Appendix contains
some additional discussion on how the choice of boundary
conditions affects the partitioning of the change in potential
energy into changes in strain energy and gravitational po-
tential energy. Therefore, in the energy balance we settle for
the total change in potential energy, which is the sum of
these two energies.
Dynamic Rupture in a Homogeneous Half-Space
We first attempt to generate realistic ruptures in a ho-
mogeneous half-space. We use some of the basic features of
the ruptures observed in nature to judge the behavior of the
simulated ruptures. Heaton (1990) and Somerville
et al.
(1999) examined the rupture behavior of several earthquakes
and found no systematic variations in the slip distributions
with depth. In other words, the slip distributions could be
approximated to first order as the changes in strain on the
fault that are independent of depth. In a homogeneous half-
space this corresponds to stress drops that are uniform with
depth. For earthquakes in a homogeneous half-space with a
uniform stress drop, we expect tapering in the slip along the
buried edges of the fault (e.g., see Parsons
et al.,
1988).
Although difficult to resolve, Heaton (1990) and Somerville
et al.
(1999) did not find any clear variations in the duration
of slip with depth. Consequently, we prefer numerical mod-
els where slip rates do not change dramatically with depth.
We shall also assume that the nominal tectonic tractions may
be derived from the application of relatively uniform strain
fields, which is consistent with the absence of any depth
dependence found in the distribution of slip from source
inversions.
Friction Model Parameters
We wish to create a relatively uniform slip distribution
with depth in a domain where the material properties do not
change with depth and the effective fault normal stresses
increase linearly with depth as a result of the overburden
pressure. For relatively homogeneous slip, the stress drop on
the fault will generally vary proportionally with the shear
modulus and the slip,
σμ
=
CD
1
,
(8)
where
C
1
is a constant that depends on the rupture dimen-
sions. The stress drop also equals the difference between the
initial shear stress,
r
0
, and the final shear stress,
r
1
.We
assume that the initial shear stress comes from a uniform
strain field, which gives
σεμ
00
=
.
(9)
For the final shear stress, we use the minimum sliding shear
stress,
σμσ
1
=−
min
.
n
(10)
Combining these equations, substituting in the expression
for the shear modulus (
l

qb
2
, where
q
denotes the mass
density and
b
denotes the shear-wave speed), and solving
for the minimum coefficient of friction yields
μ
ερβ
σ
min
()
.
=
CD
n
10
2
(11)
We now consider two end cases for the effective normal
stresses in our homogeneous half-space: uniform effective
normal stresses and lithostatic effective normal stresses.
With uniform effective normal stresses, for uniform slip we
require a uniform minimum coefficient of friction. This set
of parameters, although not physically realistic, is typically
used in dynamic rupture simulations (Day, 1982b; Olsen
et
al.,
1997; Madariaga
et al.,
1998; Harris and Day, 1999;
Oglesby
et al.,
2000b). On the other hand, if the effective
normal stress increases linearly with depth, then for uniform
slip we need a minimum coefficient of friction that varies
inversely with depth. If we try to use a uniform minimum
coefficient of friction when the normal stress increases with
depth, then the slip will increase rapidly with depth (Aa-
gaard, 1999). Such behavior conflicts with the inferred slip
1772
B. T. Aagaard, T. H. Heaton, and J. F. Hall
x
y
S1
S2
N
o
rth
H
32.00km
30
.0
0
km
1
00
.0
0k
m
60
.0
0k
m
40.00km
15.00km
20.00km
z
Figure
5.
Domain geometry for the strike-slip fault. The label H denotes the hypo-
center location.
distributions from source inversions that exhibit no clear
trends with depth (Heaton, 1990; Somerville
et al.,
1999).
The procedure just described does not yield information
regarding the maximum coefficient of friction that we as-
sociate with the failure stress and fracture energy. We choose
to vary the maximum coefficient of friction with the inverse
of the depth as well because this leads to uniform relative
changes in the coefficient of friction during sliding. This
corresponds to relatively uniform failure stresses, sliding
stresses, and fracture energies with depth in the homoge-
neous half-space.
It is important to note that these constraints on the co-
efficient of friction arise from our choice to follow the con-
ventional formulation of the friction stress (friction stress
equals the product of the coefficient of friction and the nor-
mal stress). Alternatively, a more physically meaningful
approach might be to assume that the friction force does
not depend on the normal stress. For example, the slip-
weakening friction model can be reformulated to yield the
friction stress as a function of slip with low (or zero) sliding
friction replacing the minimum coefficient of friction and a
drop in friction stress over a given slip distance acting as a
proxy for the fracture energy.
Application to Strike-Slip Fault
We now demonstrate how this parameterization works
for the case of dynamic rupture on a strike-slip fault in a
homogeneous half-space. The 60-km-long and 15-km-wide
strike-slip fault lies in a domain 100 km long, 40 km wide,
and 32 km deep as shown in Figure 5. The homogeneous
half-space has a mass density of 2450 kg/m
3
, a shear-wave
speed of 3.30 km/sec, and a dilatational wave speed of 5.70
km/sec. The finite-element model contains 3.0 million ele-
ments and 1.8 million degrees of freedom.
We assume that the effective normal stresses equal the
lithostatic pressures. In order to select the initial shear trac-
tions on the fault, we begin by choosing a maximum dy-
namic stress drop of 2.0 MPa. Recall that because the co-
efficient of friction returns to its maximum value upon
termination of sliding, the dynamic stress drop exceeds the
static stress drop. We assume that the earthquake does not
completely relieve the initial stresses and apply uniform ini-
tial shear tractions of 4.0 MPa that are tapered at the buried
edges of the fault to smother the rupture as it approaches
these edges. Figure 6 gives the distributions of the normal
and shear tractions on the fault surface.
We also need to determine a value for the failure stress
in order to specify the maximum value for the coefficient of
friction. We expect the initial stresses to lie somewhere be-
tween the minimum sliding shear stresses and the shear
stresses at failure. A small distance from failure (the differ-
ence between the failure stress and the initial shear stress)
implies that the fault is close to failure and that the rupture
will propagate very fast. However, note that the choice of
failure shear stress depends on the discretization size and
that we are actually changing fracture energies. At the other
extreme, a large distance from failure inhibits propagation
of the rupture. We avoid the extreme cases and select the
distance from failure to match the maximum dynamic stress
drop. As a result, the initial shear stresses of 4.0 MPa lie
halfway between the minimum sliding shear stresses of 2.0
MPa and the failure stresses of 6.0 MPa.
Using the normal stresses, we select values for the fric-
tion model parameters to yield these choices for the stress
drop and the distance from failure. We use the slip-weaken-
ing friction model with the two parameters,
l
max
and
l
min
,
decreasing with the inverse of depth as given by
Dynamic Earthquake Ruptures in the Presence of Lithostatic Normal Stresses
1773
15
10
5
0
Dist. down Dip (km)
Shear Traction
Traction (MPa)
0.0
1.5
3.0
4.5
6.0
0
10
20
30
40
50
60
15
10
5
0
Dist. along Strike (km)
Dist. down Dip (km)
Normal Pressure
Pressure (MPa)
0
100
200
300
400
Figure
6.
Initial shear tractions and normal pressures on the fault surface of the
homogeneous half-space. The shear tractions are derived from a uniform strain field,
and the normal pressures increase with depth because of gravity.
μ
μ
max
min
.
.
.
=
>−
<−
=
>−
1 00
250
250
250
0 333
250
83 3
z
z
z
z
z
m
m
m
m
m
z
z
D
<−
=
250
0 150
0
m
m
..
(12)
It seems unreasonable to let the coefficient of friction ap-
proach infinity at the surface, so we clip its value above a
depth of 250 m. The slip-weakening distance of 0.150 m
corresponds to a fracture energy (3

10
5
J/m
2
) that yields
reasonable rupture speeds.
The snapshots of slip rate in Figure 7 illustrate several
important features of the rupture. The rupture expands as an
ellipse with a faster rupture speed in the direction of slip
compared to the direction perpendicular to slip. In the di-
rection of slip, the rupture displays mode-II crack behavior
(shearing), and in the direction perpendicular to slip, the rup-
ture displays mode-III crack behavior (tearing). In the ab-
sence of fracture energy, mode-II cracks propagate at the
Rayleigh-wave speed, and mode-III cracks propagate at
the shear-wave speed (Freund, 1990). However, for reason-
able rupture speeds (fracture energies), the rupture speed in
the direction of slip (mode-II) tends to exceed the speed in
the direction perpendicular to slip (mode-III), because of the
asymmetry in the shear-wave radiation pattern (Andrews,
1976b; Day, 1982b; Madariaga
et al.,
1998). We observe
precisely this type of behavior as illustrated by the solid
ellipses in the figure that identify the leading edge of the
rupture; in the direction parallel to the slip we observe a
rupture speed of 2.2 km/sec (compared with a shear-wave
speed of 3.3 km/sec), and in the direction perpendicular to
the slip we observe a rupture speed of 1.8 km/sec.
When the rupture reflects off the free surface, we fold
over the solid ellipse at 5.5 sec to coincide with the reflected
portion of the rupture. We also add a dashed ellipse at the
leading edge of another portion that propagates along the
free surface at a speed of 4.4 km/sec. This speed lies between
the shear-wave speed of 3.3 km/sec and the dilatational-
wave speed of 5.7 km/sec. The rupture at the free surface
begins to separate at 8.5 sec, and at 10 sec two distinct slip
events occur at the free surface. We observe substantially
larger slip rates where the two portions of the rupture con-
structively interfere (identified by the intersection of the
solid and dashed ellipses). In general, the portion of the rup-
ture traveling faster than the shear-wave speed is associated
with larger slip rates than the portion traveling below the
shear-wave speed. Furthermore, as the two portions interact,
the speed of the portion that propagates below the shear-
wave speed increases to around 3.0 km/sec. The solid ellipse
that we overlay on the slip rate snapshots reflects this change
in rupture speed.
This complex rupture yields an average slip of 1.0 m
and creates the smooth distribution of final slip shown in
Figure 8, which is what we expect based on the uniform
maximum dynamic stress drop. In other words, beyond the
tapering along the edges of the fault, the slip distribution does
not exhibit any clear trends with depth, so it is compatible
with those from source inversions. In contrast with the final
slip, the peak slip rate reflects the complex nature of the rup-
1774
B. T. Aagaard, T. H. Heaton, and J. F. Hall
15
10
5
0
Time = 1.0 sec
15
10
5
0
Time = 2.5 sec
15
10
5
0
Dist. down Dip (km)
Time = 4.0 sec
15
10
5
0
Time = 5.5 sec
0
10
20
30
40
50
60
15
10
5
0
Dist. along Strike (km)
Time = 7.0 sec
Time = 8.5 sec
Time = 10.0 sec
Time = 11.5 sec
Time = 13.0 sec
0
10
20
30
40
50
60
Dist. along Strike (km)
Time = 14.5 sec
Slip Rate (m/sec)
0.0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1.0
Figure
7.
Snapshots of slip rate on the fault surface in the homogeneous half-space
with the slip-weakening friction model. The solid and dashed ellipses indicate the
leading edges of the ruptures propagating slightly below the shear-wave speed, and
between the shear- and dilatational-wave speeds, respectively.
ture. The path of constructive interference between the two
portions of the rupture is clearly visible, with slip rates
roughly 0.2 m/sec greater than in the surrounding regions. We
also see large slip rates near the top of the fault at the north
end, which come from the faster portion of the rupture.
The faster portion of the rupture propagates at a speed
between the shear-wave speed and the dilatational-wave
speed, whereas the slower portion propagates at a speed be-
low the shear-wave speed. Both Burridge
et al.
(1979) and
Freund (1979) found steady-state solutions for propagation
at speeds both slower than the Rayleigh-wave speed and
between the shear-wave speed and the dilatational-wave
speed. They concluded that stable propagation can occur for
mode-II cracks with speeds of
2
times the shear-wave
speed. Furthermore, Rosakis
et al.
(1999) observed cracks
propagating at
2
times the shear-wave speed in a brittle
polyester resin under far-field loading. In our simulation, the
faster portion of the rupture propagates at approximately 1.3
times the shear-wave speed or within 6% of
2
times the
shear-wave speed.
The super-shear rupture speeds in our simulation arise
from the large slip rates along the ground surface. When the
rupture hits the free surface, it encounters a reduced resis-
tance to slip. With a hypocenter several kilometers below
the surface, a portion of the rupture front several kilometers
long hits the ground surface nearly simultaneously and cre-
ates a high apparent velocity along the ground surface. Con-
sequently, the slip rates increase and the rupture sustains a
super-shear propagation speed along the surface. On a long
fault (as in our strike-slip fault above) this leads to bifurca-
tion of the rupture into a portion propagating at approxi-
mately
2
times the shear-wave speed and a portion prop-
agating a little below the shear-wave speed.
Whereas the portion propagating faster than the shear-
wave speed appears feasible, source inversions of earth-
quakes indicate that the ruptures generally propagate at
Dynamic Earthquake Ruptures in the Presence of Lithostatic Normal Stresses
1775
15
10
5
0
Dist. down Dip (km)
Final Slip
Slip (m)
0.0
0.5
1.0
1.5
2.0
0
10
20
30
40
50
60
15
10
5
0
Dist. along Strike (km)
Dist. down Dip (km)
Peak Slip Rate
Slip Rate (m/sec)
0.00
0.25
0.50
0.75
1.00
Figure
8.
Distributions of final slip and peak slip rate on the fault surface in the
homogeneous half-space. The bifurcation of the rupture generates complexity in the
distribution of peak slip rate but has little effect on the distribution of final slip.
speeds below the shear-wave speed (Heaton, 1990; Somer-
ville
et al.,
1999), although some evidence exists for super-
shear rupture speeds over small portions of fault ruptures
(Archuleta, 1984; Hernandez
et al.,
1999). This discrepancy
in the speed of propagation may be explained by two short-
comings of the numerical simulations of ruptures in a dis-
cretized, homogeneous half-space. A finer discretization size
allows a larger failure stress for the same fracture energy.
This tends to localize the stress concentration and to inhibit
the transition to super-shear rupture speeds. Higher fracture
energies would also impede the transition to supershear rup-
ture speeds. Additionally, the Earth includes variations in
the material properties with generally softer material at shal-
low depths. This reduces the initial stress on the fault for a
given amount of strain and reduces the rupture speed as it
propagates toward the surface. Hence, the use of too small
a fracture energy for the homogeneous material properties
may be responsible for the discrepancy in the rupture speeds
between the simulation and the real earthquakes.
Dynamic Rupture in a Layered Half-Space
Having created a dynamic rupture on a strike-slip fault
in a homogeneous half-space that exhibits behavior gener-
ally compatible with source inversions, we now study dy-
namic rupture on faults in a layered half-space. Additionally,
we attempt to match the general characteristics of the dy-
namic ruptures with those of the prescribed ruptures from
our previous study of near-source ground motions (Aagaard,
1999; Aagaard
et al.,
2001).
We consider both a vertical strike-slip fault and a
shallow-dipping thrust fault. The geometry of the strike-slip
domain matches that of the strike-slip fault in the homoge-
neous half-space discussed earlier and shown in Figure 5.
The geometry of the thrust fault resembles the geometry of
the Elysian Park fault underneath Los Angeles (Hall
et al.,
1995). The 28-km-long and 18-km-wide thrust fault dips 23

to the north. Figure 9 shows the geometry of the thrust fault
with the top of the fault sitting 8.0 km below the ground
surface. The domains for both faults contain the same vari-
ation of the material properties with depth. Figure 10 gives
the mass density, shear-wave speed, and dilatational-wave
speed as a function of depth in the layered half-space. We
partition the finite-element meshes among 16 processors us-
ing the METIS library (Karypis
et al.,
1999) from the Uni-
versity of Minnesota. The strike-slip fault simulations re-
quired 5.6 hr and the thrust fault simulations required 2.6 hr
on the Hewlett-Packard X-Class supercomputer at the Center
for Advanced Computing Research at Caltech.
Friction Model Parameters
We return to equation (11),
μ
ερβ
σ
min
()
,
=
CD
n
10
2
(13)
with the objective of creating slip distributions that do not
vary systematically with depth. In contrast with the homo-
geneous half-space, the mass density and shear-wave speed
follow a complex variation with depth in the layered half-
1776
B. T. Aagaard, T. H. Heaton, and J. F. Hall
S1
North
60.00km
60.00km
23
°
S2
8.00km
18.00km
28.00km
25.00km
24.00km
z
y
x
H
Figure
9.
Domain geometry for the thrust fault for the case where the top of the
fault lies 8.0 km below the ground surface. The label H denotes the hypocenter location.
0
1
2
3
4
5
6
7
30
20
10
0
Velocity (km/s), Density (g/cm
3
)
Depth (km)
α
β
ρ
Figure
10.
Density (
q
), shear-wave speed (
b
), and
dilatational-wave speed (

) as a function of depth in
the layered half-space.
space. Above a depth of 6.0 km, the values increase with
depth piecewise linearly, whereas below a depth of 6.0 km
the material properties are relatively uniform. Making the
coefficient of friction proportional to either the ratio of the
square root of the shear modulus to the depth or the ratio of
the shear-wave speed to the depth provides a reasonable
match to the desired variation of
l
min
given by equation (13),
with the material properties in the layered half-space. We
choose to vary the minimum coefficient of friction with the
ratio of the square root of the shear modulus to the depth.
As noted earlier in our discussion of the friction model pa-
rameters for the homogeneous half-space, formulating the
friction stress from the product of the coefficient of friction
and the normal stress in the presence of lithostatic normal
pressures requires variations in the coefficient of friction
over the depth of the fault.
Simulation Nomenclature
We name each scenario based on the choice of simu-
lation parameters. Table 2 displays the correspondence be-
tween the letters and numbers of the scenario names and the
simulation parameters.
Application to Strike-Slip Fault
We follow the same procedures that we used for the
dynamic rupture simulations in the homogeneous half-space
to determine the initial shear and normal tractions on the
fault surface in the layered half-space. For each scenario, we
assume that the pore pressures are negligible and apply litho-
static effective normal tractions. We aim for an average slip
of 2.0 m. The average final stress drop on a rectangular,
vertical strike-slip fault in a homogeneous Poissonian half-
space approximately follows
σμ
=
CDw
(/),
where
C
Cwl
Clwwlw
Clw
C
D
D
D
D
=
+−<<
−>
091
2
09
2
16
.[ (/ )]
.
.
for surface
e rupture
for deeply buried faults,
21
.
(14)
where
l
and
w
denote the length and width of the fault (Hea-
ton
et al.,
1986). Below a depth of 6.0 km, the material
properties on the fault surface are nearly uniform, so we use
the shear modulus from a depth of 6.0 km in equation (14).
Applying this equation with an average slip of 2.0 m and
our fault dimensions yields an average stress drop of 2.5
MPa. The recovery of the coefficient of friction upon ter-
mination of sliding means that the maximum dynamic stress
drop will exceed the average final stress drop. Consequently,
Dynamic Earthquake Ruptures in the Presence of Lithostatic Normal Stresses
1777
Table 2
Description of Letters and Numbers Used to Compose the Scenario Names*
Fault Type
Fault Depth
/ Friction Model
/ Heterogeneity
SS (strike slip)
0 (0 km)
SW (slip weakening)
U (uniform)
TH (thrust)
8 (8 km)
SRW (slip and rate weakening)
HS (heterogeneous stress/strain)
P (prescribed rupture)
HF (heterogeneous friction)
*For example, scenario SS0/SW/U refers to a strike-slip scenario with the top of the fault at a depth of 0
km, a slip-weakening friction model, uniform initial shear tractions, and uniform coefficients in the expression
for the coefficient of friction
on the basis of a test simulation with a homogeneous half-
space, we impose a maximum dynamic stress drop of 4.5
MPa below a depth of 6.0 km, where the material properties
are relatively uniform.
We assume that some residual shear stresses remain on
the fault after the earthquake and that the initial tectonic
shear tractions come from a relatively uniform strain field,
so we select tectonic shear strains that produce shear trac-
tions of approximately 6.0 MPa at depths where the material
properties are nearly uniform. In order to prevent the effec-
tive normal stress on the fault surface from vanishing at the
ground surface, we also apply uniform axial tectonic strains
in the east–west direction (parallel to the
x
axis). The tectonic
strain field (denoted by the superscript
t
) is given by
ε
εενε
ε
εε
xx
t
yy
t
zz
t
xx
t
xy
t
yz
t
xz
t
=−×
==−
==
293 10
110 10
0
5
4
.
.
,
(15)
where
m
denotes Poisson’s ratio. We superimpose this strain
field on the strains associated with the lithostatic tractions
generated by gravity (denoted by the superscript
g
),
εεε
λμ
ρ
εεε
xx
g
yy
g
zz
g
z
xy
g
yz
g
xz
g
sg s
===
+
===
1
32
0
0
()
,
d
(16)
to generate the shear and normal tractions that we apply on
the fault surface (Fig. 11).
Having chosen to vary the minimum coefficient of fric-
tion with the square root of the shear modulus and the in-
verse of the depth, we apply the same functional form to the
maximum coefficient of friction. This creates a uniform rela-
tive change in the coefficient of friction with depth that we
also used in the homogeneous half-space. With the varia-
tions of the material properties with depth, we cannot match
the distance from failure to the maximum dynamic stress
drop over the entire depth of the fault as we did in the ho-
mogeneous half-space. As a result, we choose to match (in
an average sense) the distance from failure with the maxi-
mum dynamic stress drop over the depth range of 6.0–15.0
km, where the material properties remain relatively uniform.
The functional forms of the parameters in the slip-weakening
friction model are given by
μ
μ
max
/
..
..
=
>−
−×
<−
0 164
1 0
302 10
10
3
2
12
z
z
z
km
msec
kg
km
3
=
>−
−×
<−
μ
μ
min
/
..
..
0 0235
1 0
431 10
10
4
12
z
z
z
km
msec
kg
k
32
m
m
m
=
D
0
0 338
..
(17)
We clip the values above a depth of 1.0 km to prevent them
from approaching infinity at the ground surface. Note that
l
max
and
l
min
denote constants in the friction model, whereas
l
denotes rigidity. Figure 12 shows the maximum and min-
imum coefficients of friction as a function of distance down-
dip on the fault.
Figure 12 also displays the initial shear stress through
the center of the asperity used to start the rupture, the shear
stress at failure, and the minimum sliding shear stress over
the depth of the fault. The maximum dynamic stress drop
(difference between the initial and minimum sliding shear
stresses) closely follows the variation of the shear modulus
that increases piecewise linearly down to a depth of 6.0 km
and is nearly uniform below 6.0 km; this means that the
change in shear strain, and hence of slip, will tend to be
relatively uniform. The fracture energy per unit area (Fig.
12) follows the variation of the failure and the minimum
sliding shear stresses with depth. Consequently, the fracture
energy remains nearly independent of depth below a depth
of 6.0 km, but becomes progressively smaller above a depth
of 6.0 km. Using the nomenclature given in Table 2, we shall
refer to this scenario as SS0/SW/U.
The rupture begins propagating at about 2.5 km/sec in
the direction parallel to slip (mode-II), and after hitting the
ground surface it maintains a speed of 3.0 km/sec (91%
of the local shear-wave speed) in the direction parallel to
the slip at a depth of 6.0 km. We attribute the change in the
rupture speed to the increase in the peak slip rates as the
rupture approaches the ground surface and encounters a re-
duction in the stiffness and fracture energy. The rupture re-
flects off the ground surface, but this additional slip soon
1778
B. T. Aagaard, T. H. Heaton, and J. F. Hall
15
10
5
0
Dist. down Dip (km)
Shear Traction
Traction (MPa)
0.0
3.0
6.0
9.0
12.0
0
10
20
30
40
50
60
15
10
5
0
Dist. along Strike (km)
Dist. down Dip (km)
Normal Pressure
Pressure (MPa)
0
100
200
300
400
Figure
11.
Initial shear tractions and normal pressures on the fault for scenario
SS0/SW/U. The tractions result from the superposition of the tectonic strains and the
strains caused by gravity (equations 15 and 16).
0.00
0.05
0.10
0.15
0.20
15
12
9
6
3
0
Coefficient of Friction
Dist. down Dip (km)
μ
max
μ
min
0
3
6
9
12
Shear Stress (MPa)
initial
failure
sliding
0.0
0.5
1.0
1.5
2.0
Fracture Energy (MJ/m
2
)
Figure
12.
Maximum and minimum coefficients of friction (left), initial, failure,
and minimum sliding shear stresses through the center of the asperity (center), and
fracture energy per unit area (right) as functions of distance downdip on the fault for
scenario SS0/SW/U. The fracture energy per unit area follows the variation of the
failure and minimum sliding shear stresses with depth. Whereas we can increase
l
max
significantly without altering the rupture behavior by decreasing the discretization size,
increasing
l
min
results in large temperature changes on the fault (see Fig. 15).
disappears. In contrast with the homogeneous half-space
simulation, the rupture does not bifurcate and propagates
slower than the local shear-wave speed. The average slip of
1.9 m nearly matches the target value of 2.0 m and corre-
sponds to a moment magnitude of 6.9. We compute a final
stress drop (averaged over the fault surface) of 1.4 MPa,
which falls short of the 2.5 MPa final, uniform stress drop
predicted by equation (14) for a homogeneous half-space.
The presence of the softer material in the top 6.0 km of the
domain reduces the average stress drop compared with that
of a homogeneous half-space with the same average slip.
Figure 13 shows the distributions of the final slip and
peak slip rate on the fault surface. The region where the final
slip exceeds 3.0 m coincides with the locations subjected to
the additional slip associated with the reflection of the rup-
ture off the free surface. We find that the peak slip rates near
Dynamic Earthquake Ruptures in the Presence of Lithostatic Normal Stresses
1779
15
10
5
0
Dist. down Dip (km)
Final Slip
Slip (m)
0.0
1.0
2.0
3.0
4.0
0
10
20
30
40
50
60
15
10
5
0
Dist. along Strike (km)
Dist. down Dip (km)
Peak Slip Rate
Slip Rate (m/sec)
0.0
0.5
1.0
1.5
2.0
Figure
13.
Distributions of final slip and peak slip rate on the fault for scenario
SS0/SW/U. Both the final slip and the peak slip rate are, to first order, independent of
depth after removing the tapering along the buried edges of the fault.
the surface are about 0.5 m/sec greater than the peak slip
rates at depth. The slight tendency for the peak slip rates to
increase as the rupture propagates (the peak slip rates in-
crease by roughly 0.25 m/sec over a distance of 25 km)
causes this region of larger slip rates at the surface to pro-
gressively increase in size. In contrast with Mikumo (1992),
who also used normal stresses equal to the overburden pres-
sure, the slip distributions here do not exhibit a strong vari-
ation with depth. This difference stems from the fact that we
chose the friction parameters such that the dynamic stress
drop follows the change in shear modulus in order to pro-
duce relatively uniform changes in strain (slip). Mikumo
(1992) decreased the dynamic stress drop with depth follow-
ing the transition from brittle to ductile behavior of materi-
als, and, as a result, the slip decreased dramatically with
depth.
The peak horizontal particle displacements and veloci-
ties on the ground surface in Figure 14 clearly illustrate the
effect of the rupture directivity on the ground motions. Both
the peak horizontal displacements and the velocities increase
along the strike of the fault away from the epicenter until
the end of the fault where they begin to decay. The peak
displacements exceed 1.0 m over an area of approximately
1200 km
2
with a maximum value of 3.0 m. Likewise, the
peak velocities exceed 1.0 m/sec over an area of approxi-
mately 550 km
2
with a maximum value of 3.5 m/sec. Al-
though the peak displacements decay away from the fault at
a slower rate than the peak velocities, the most severe motion
is confined to a narrow region along the fault. These near-
source ground motions display the same principal features
as those of Olsen and Archuleta (1996) and Olsen
et al.
(1997) for similar-sized events on strike-slip faults.
We now evaluate the level of the shear stresses on the
fault during sliding. Recall that in most instances adding a
constant value to the initial, failure, and sliding shear stresses
does not significantly change the rupture behavior (Guatteri
and Spudich, 1998). This means that the observed rupture
behavior and ground motions do not constrain the absolute
levels of these stresses. However, using equation (6) we can
compute the change in temperature at each point on the fault.
Observations of exposed fault surfaces indicate that slip in
an earthquake likely occurs across a region of less than a
few millimeters with little melting of the rocks (Chester and
Chester, 1998). This implies that the change in temperature
on the fault during sliding remains below the level that
would cause melting. This is a variation of the classic heat-
flow problem on the San Andreas fault (Brune
et al.,
1969;
Lachenbruch, 1980). In our simple analysis, we shall assume
that the initial temperature at each point on the fault is small
compared with the melting temperature and that changes in
temperature of the order of 1000 K cause melting. Thus, we
want our sliding stresses to produce temperature changes
less than 1000 K.
We estimate the temperature change on the fault during
sliding using equation (6) by following a procedure similar
to that of McKenzie and Brune (1972), Richards (1976), and
Kanamori
et al.
(1998). We begin by assuming a heat ca-
pacity per unit mass of 1000 J/kg K. If the slip occurs across
1780
B. T. Aagaard, T. H. Heaton, and J. F. Hall
0
10
20
30
40
50
60
15
10
5
0
Dist. along Strike (km)
Dist. down Dip (km)
Change in Temp. (
o
K)
0
100
200
300
400
Figure
15.
Final change in temperature on the fault for scenario SS0/SW/U. By
using very low values of sliding friction (see Fig. 12), the temperature changes are
compatible with the lack of melting observed in exposed faults.
20

10
0
10
20
East–West (km)
Peak Horizontal Displacement
Ma
g
nitude
(
m
)
0.0
0.5
1.0
1.5
2.0
2.5
3.0
20

10
0
10
20
50
40
30
20
10
0
10
20
30
40
50
East–West (km)
South–North (km)
Peak Horizontal Velocity
Ma
g
nitude
(
m/sec
)
0.0
0.5
1.0
1.5
2.0
2.5
3.0


Figure
14.
Peak horizontal particle dis-
placements and velocities on the ground sur-
face for scenario SS0/SW/U. The line indicates
the projection of the fault plane on to the
ground surface, and the hollow circle identifies
the epicenter. The amplitude of the ground mo-
tions increases along the fault away from the
epicenter because of rupture directivity.
an infinitesimally thin zone, then the heat is confined to the
thermal penetration depth given by
dk
d
=
τ
,
where
k
is the
thermal diffusivity and
s
d
the timescale of the slip. Assuming
that
k

1.35

10

6
m/sec
2
and choosing
s
d

5 sec gives
d

2.6 mm. If the slip is distributed over a zone of finite
width, then
d
would be larger. Consequently, we moderate
the value of
d
given by the thermal penetration depth with
an infinitesimally thin slip zone and choose a value of
d

5.0 mm. As shown in Figure 15, at most locations the tem-
perature increases by 200–300 K. We observe smaller
changes in temperature near the top of the fault because the
sliding stresses are smaller there. Below a depth of 6.0 km
the sliding stresses vary little; therefore, below that depth the
change in temperature closely resembles the distribution of
final slip. Thus, our sliding stresses seem consistent with the
lack of melting observed in fault-zone materials. However,
if the minimum coefficient of friction was increased to levels
reported from laboratory measurements (
l
min

0.6 [Dieter-
ich, 1992]), then the temperature changes at middepth on
the fault would be 64 times larger than those shown in
Figure 15.
On the basis of this scenario, we find that uniform tec-
tonic strains and a friction model with parameters that vary
proportionally to the square root of the shear modulus and
inversely with depth produce a realistic rupture: we observe
rupture speeds and slip distributions that are compatible with
source inversions, the ground motions exhibit the directivity
that we expect, and the estimated changes in temperature on
the fault surface remain consistent with the limited amount
of glassy material observed in fault zones. Moreover, chang-
ing the dependence of the coefficient of friction from the
square root of the shear modulus to the shear-wave speed
also yields nearly identical behavior (Aagaard, 1999). Dy-
namic ruptures with uniform effective normal stresses dis-
Dynamic Earthquake Ruptures in the Presence of Lithostatic Normal Stresses
1781
play similar behavior (e.g., Day 1982b; Mikumo and Mi-
yatake, 1993; Madariaga
et al.,
1998) but would require
extremely high pore pressures at depth to be physically
meaningful.
Application to Thrust Fault
We now turn our attention to dynamic ruptures on the
thrust fault in the layered half-space (Fig. 9). On dipping
faults, the tectonic stresses generate both shear and normal
stresses on the fault surface. With pore pressures at or below
the hydrostatic pressure, gravity creates effective normal
stresses that increase with depth and far exceed the contri-
butions of the tectonic stresses’ to the normal stresses. Con-
sequently, changing the dip angle of the fault while keeping
the tectonic stresses constant causes almost no variations in
the effective normal stresses. This means we can follow the
same procedure that we used for the strike-slip fault for dy-
namic failure on the thrust fault. This would not be the case
if we chose to neglect the effect of gravity on the normal
stresses and used uniform effective normal stresses.
We apply uniform horizontal axial strains and shear
strains to generate the shear tractions on the fault. We align
the shear tractions with the desired slip direction, which has
a rake angle of 105

from the strike, and aim for an average
slip of 1.0 m. For inclined faults the average stress drop
remains proportional to the product of the shear modulus
and average slip, but the proportionality constant depends
on the depth and dip angle of the fault. Consequently, for
inclined faults we do not have a simple expression for the
average stress drop as a function of the shear modulus and
the average slip that we have for strike-slip faults (equation
14) (Parsons
et al.,
1988).
With the top of the fault buried 8.0 km below the ground
surface, the material properties exhibit little change over the
depth of the fault. As a result, uniform tectonic strains create
nearly uniform shear and normal tractions on the fault. The
shallow dip of the fault causes the tectonic strains to produce
much smaller normal tractions than the normal tractions
from gravity. We do not change the functional form of the
slip-weakening friction model from the one used in the
strike-slip case; the parameters in the friction model continue
to depend on the ratio of the square root of the shear modulus
to the depth. We do change the coefficients slightly to create
the desired maximum dynamic stress drop and shear stresses
at failure.
We use homogeneous initial tectonic strains to generate
nominal shear tractions of 6.0 MPa on the fault surface and
superimpose these tractions on the lithostatic tractions gen-
erated by gravity. We again denote the tectonic strains,
ε
εενε
ε
εε
yy
t
xx
t
zz
t
yy
t
xy
t
yz
t
xz
t
==−
==
236 10
727 10
0
4
5
.
.
,
(18)
with the superscript
t
and the strains caused by gravity,
εεε
λμ
ρ
εεε
xx
g
yy
g
zz
g
z
xy
g
yz
g
xz
g
sg s
===
+
===
1
32
0
0
()
,
d
(19)
with the superscript
g.
We select nominal minimum sliding
stresses of 1.5 MPa, a nominal maximum dynamic stress
drop of 4.5 MPa, and nominal shear stresses at failure of
10.5 MPa, which yields
μ
μ
max
/
..
.
sec
.
=
>−
−×
<−
0 162
1 0
297 10
10
3
32
12
z
z
z
km
m
kg
km
=
>−
−×
<−
μ
μ
min
/
..
.
sec
.
0 0231
1 0
424 10
10
4
32
12
z
z
z
km
m
kg
km
m
m
=
D
0
0 338
.
(20)
for the parameters in the slip-weakening friction model.
These match the stresses on the strike-slip fault at similar
depths. We will refer to this scenario as TH8/SW/U.
We start the rupture using a shear stress asperity with a
radius of 1.8 km located along the north-south centerline of
the fault at a depth of 13.5 or 4.0 km updip from the bottom
of the fault (Fig. 9, hypocenter H). The taper in the shear
tractions on all four edges smothers the rupture as it ap-
proaches the edges of the fault. Figure 16 displays the initial
shear and normal tractions applied to the fault surface.
The rupture begins slowly in response to the placement
of the asperity close to the edge of the fault. As the rupture
begins to propagate, the rupture front conforms to the fa-
miliar elliptic shape with the fastest rupture speed in the
direction of slip, which has a rake angle of 105

. The peak
slip rates remain relatively low, and the rupture propagates
in the direction of slip at a speed of only 2.2 km/sec or about
67% of the local shear-wave speed. The elliptic shape of the
rupture front causes the leading edge of the rupture to reach
the top of the fault several seconds before the rupture reaches
the lateral edges; this gives the rupture a bilateral appear-
ance.
The distribution of final slip displayed in Figure 17 dis-
plays no clear trends with depth and resembles the final slip
of a statically applied uniform stress drop. The average slip
of 1.2 m agrees reasonably well with the target value of 1.0
m. The reflection of the dilatational wave off the ground
surface generates a shear wave that propagates back down
toward the fault. As the wave passes through the fault, the
dynamic shear stresses cause additional sliding on the fault
and the peak slip of 2.3 m near the hypocenter. The slip rates
associated with this additional slip near the hypocenter ex-
ceed those in the same region for the first slip event, which
reflect the slow initiation of the rupture. If we neglect the