of 18
Constraining fault constitutive behavior with slip
and stress heterogeneity
B. T. Aagaard
1
and T. H. Heaton
2
Received 10 October 2006; revised 5 November 2007; accepted 16 January 2008; published 8 April 2008.
[
1
]
We study how enforcing self-consistency in the statistical properties of the preshear
and postshear stress on a fault can be used to constrain fault constitutive behavior beyond
that required to produce a desired spatial and temporal evolution of slip in a single
event. We explore features of rupture dynamics that (1) lead to slip heterogeneity in
earthquake ruptures and (2) maintain these conditions following rupture, so that the stress
field is compatible with the generation of aftershocks and facilitates heterogeneous slip in
subsequent events. Our three-dimensional finite element simulations of magnitude 7
events on a vertical, planar strike-slip fault show that the conditions that lead to slip
heterogeneity remain in place after large events when the dynamic stress drop (initial shear
stress) and breakdown work (fracture energy) are spatially heterogeneous. In these
models the breakdown work is on the order of MJ/m
2
, which is comparable to the radiated
energy. These conditions producing slip heterogeneity also tend to produce narrower slip
pulses independent of a slip rate dependence in the fault constitutive model. An alternative
mechanism for generating these confined slip pulses appears to be fault constitutive
models that have a stronger rate dependence, which also makes them difficult to
implement in numerical models. We hypothesize that self-consistent ruptures could also be
produced by very narrow slip pulses propagating in a self-sustaining heterogeneous stress
field with breakdown work comparable to fracture energy estimates of kJ/M
2
.
Citation:
Aagaard, B. T., and T. H. Heaton (2008), Constraining fault constitutive behavior with slip and stress heterogeneity,
J. Geophys. Res.
,
113
, B04301, doi:10.1029/2006JB004793.
1. Introduction
[
2
] Uncovering the physics controlling earthquake rup-
ture is an inherently difficult task because ruptures are
relatively infrequent and hidden from direct observation
several kilometers below the ground surface. Nevertheless,
with current computational resources it is becoming much
more common for researchers to strive to reproduce the
dynamics of natural earthquakes. When attempting to tie
models of earthquake rupture to natural behavior we must
keep in mind that there is likely a range of model parameters
that produces similar results and a model may capture only
the most general physics of the rupture process. Under-
standing the ramifications of this second issue requires a
thorough knowledge of the numerical limitations of a model
and the range of physical processes capable of generating
seismologic and geodetic observations.
[
3
] Current numerical models of spontaneous dynamic
ruptures involve spatial and temporal discretizations of the
model space, using techniques such as finite differences
[e.g.,
Andrews
, 1976a;
Day
, 1982], finite elements [e.g.,
Mikumo and Miyatake
, 1978;
Oglesby et al.
, 1998;
Aagaard
et al.
, 2001], or boundary integral techniques [e.g.,
Lapusta
et al.
, 2000]. As a result, the range of length scales, even
with state-of-the-art supercomputers, is limited to a few
orders of magnitude. Earthquake ruptures, on the other
hand, are often associated with fractal distributions (e.g.,
Gutenberg-Richter frequency-magnitude relation and fractal
spatial distributions of faults and earthquakes) [
Turcotte and
Malamud
, 2002]. In terms of rupture behavior, one could
infer from the fractal nature of earthquake behavior that slip
occurs via slip pulses in small events in much the same
fashion as it does in large events. In other words, the length
scale of the slip pulse (i.e., the distance from the leading
edge of the rupture to the trailing edge of the rupture at any
given instant during the earthquake) could be several orders
of magnitude smaller in size than the total rupture length.
[
4
] Very narrow, confined slip pulses would also be
compatible with sudden changes in sliding friction from
high static values (coefficient of friction,
m
f

0.6) to low
dynamic values (
m
f
< 0.1) and back again. Faults would
appear to be ‘‘strong’’ in the static sense but ‘‘weak’’ in the
dynamic sense, solving the heat paradox associated with the
absence of ubiquitous heat production on faults from
earthquake ruptures [
Heaton
, 1990]. Fault constitutive
models with a strong rate dependence produce such narrow
slip pulses [e.g.,
Cochard and Madariaga
, 1994]. Unfortu-
nately, they lead to ill-posed numerical solutions, either due
to the absence of a breakdown zone (for the case of pure
rate-dependent friction) or a breakdown zone much smaller
JOURNAL OF GEOPHYSICAL RESEARCH, VOL. 113, B04301, doi:10.1029/2006JB004793, 2008
1
U.S. Geological Survey, Menlo Park, California, USA.
2
Department of Geologic and Planetary Sciences, California Institute of
Technology, Pasadena, California, USA.
Copyright 2008 by the American Geophysical Union.
0148-0227/08/2006JB004793$09.00
B04301
1of18
than what can be resolved with current discretization
schemes [
Rice
, 1993;
Madariaga and Cochard
, 1994].
For example, friction stress changes occurring over a few
milliseconds (which, for the linear slip-weakening friction
model, correspond to values of slip-weakening parameter in
the range of a few millimeters) require discretization sizes
on the order of a few meters, which is nearly 2 orders of
magnitude smaller than what most current three-dimensional
(3-D) numerical models of spontaneous rupture use. Thus,
our current numerical models do not span a sufficient range
of length scales to simultaneously capture these very short
length scales while still capturing the large-scale 3-D effects
of fault slip and rupture propagation.
[
5
] One possible mechanism for working around this
difficult issue of capturing the vast range of length scales
involved in dynamic ruptures is to assume that we can
encapsulate small length-scale behavior within the fault
constitutive model, thereby permitting use of a numerical
simulation that spans only a subset of the length scales that
may be involved. We can also use such simulations to guide
our inferences in the behavior across a broader range of
scales. Keeping these limitations of the numerical models in
mind and our desire to ultimately make inferences about the
behavior across a broader range of scales, we now focus on
trying to understand the constraints on fault constitutive
models in spontaneous dynamic ruptures simulations.
[
6
] Researchers often attempt to constrain their sponta-
neous rupture models by matching strong ground motion
records or the rupture propagation imaged in kinematic
source inversions. Unfortunately, these provide only the
most basic constraints on stress drop and breakdown work
or fracture energy [see, e.g.,
Guatteri and Spudich
, 2000;
Peyrat et al.
, 2001;
Mikumo et al.
, 2003;
Zhang et al.
, 2003;
Peyrat et al.
, 2004]. In some cases [e.g.,
Guatteri and
Spudich
, 1998] these techniques have been used to con-
strain the absolute level of shear stress, but in most cases the
modeler selects the absolute level of shear stress arbitrarily.
[
7
] On the other hand, the spectral characteristics of slip
appear to be better constrained [
Somerville et al.
, 1999;
Mai
and Beroza
, 2002;
Lavalle
́e and Archuleta
, 2003;
Liu-Zeng
et al.
, 2005]. The fractal nature of the slip appears to be
associated with other scaling relationships as well, such as
the ratio of slip to rupture length and average stress drop
[
Liu-Zeng et al.
, 2005]. Several recent studies of spontane-
ous dynamic rupture have tried to exploit this to examine
the physics of the rupture, including how variations in slip,
risetime, and rupture speed affect near-source ground
motions [
Oglesby and Day
, 2002;
Guatteri et al.
, 2003]
and whether a contrast in elastic properties across a fault
affects rupture propagation in the presence of a heteroge-
neous stress field [
Andrews and Harris
, 2005]. In generating
the appropriate spatial-temporal evolution of slip, there is
often a trade-off between influencing the slip via spatial
variations in initial stress and spatial variations in the fault
constitutive properties [
Peyrat et al.
, 2004].
[
8
]
Peyrat et al.
[2004] found that they could match
strong motion records along with GPS and InSAR measure-
ments of ground deformation equally well with either (1) a
spatially heterogeneous distribution of initial shear stress
and spatially homogeneous fault constitutive parameters or
(2) a spatially homogeneous distribution of initial shear
stress and spatially heterogeneous fault constitutive param-
eters. Whereas ruptures were nearly identical from the view
of seismologic and geodetic observations (the ground
motions and static displacements they produced) the differ-
ent fault constitutive parameterizations imply the physics
controlling the ruptures were quite different. Furthermore,
the postrupture shear stress fields over the rupture area were
diametrically opposite of each other. In the case of the
initially heterogeneous shear stress field, the postrupture
shear stress field was rather homogeneous, whereas in the
case of the initially homogeneous shear stress field, the
postrupture shear stress field was quite heterogeneous. In
other words, in one case the stress field went from hetero-
geneous to homogeneous, and in the other case it went from
homogeneous to heterogeneous. One expects that natural
faults operate somewhere between these two extremes
[
Rivera and Kanamori
, 2002].
[
9
] In this study, we attempt to find sets of parameters
conducive to this self-consistent behavior associated with
faults operating between the two extremes of (1) ruptures
removing all heterogeneity from the shear stress field and
(2) ruptures introducing all heterogeneity into the shear
stress field. Heterogeneous final slip in an earthquake
corresponds to spatial variations in the static stress drop.
These spatial variations in stress drop likely develop because
both the prerupture shear stress field and the fault friction are
spatially heterogeneous. Using spontaneous dynamic rupture
simulations with different combinations of initial conditions
and fault friction parameters we try to produce ruptures with
heterogeneous slip and similar slip time histories. We then
judge whether the conditions following a rupture are consis-
tent with those preceding rupture. In this way we seek to find
additional constraints on the parameters in spontaneous
rupture simulations beyond those necessary to produce a
given spatial-temporal evolution of slip.
2. Methodology
[
10
] We study spontaneous dynamic rupture propagation
on a vertical, throughgoing strike-slip fault that is embedded
in a 75 km long, 30 km wide, and 30 km deep domain as
shown in Figure 1. Table 1 lists the principal parameters
s3
5 km
5 km
s2
15 km
s1
seismogenic zone
y
z
x
30 km
15 km
75 km
30 km
Site A
Site B
Figure 1.
Geometry of the vertical strike-slip fault and
simulation domain. The random spatial variations are tapered
outside of the seismogenic zone using equation (9), where
s
is
s
1
,
s
2
,or
s
3
. We compare velocity time histories at sites A (
x
=
0 km,
y
= 26 km) and B (
x
= 5 km,
y
=

20 km).
B04301
AAGAARD AND HEATON: FAULT FRICTION AND SLIP HETEROGENEITY
2of18
B04301
used in this study. In order to create ruptures with some of
the basic features present in natural earthquakes, we include
variation in the material properties with depth and overbur-
den pressure (hydrostatic pore pressures). The material
properties follow horizontally averaged values of
v
p
,
v
s
,
and density (given in Table 2 and Figure 2) from the
Southern California Earthquake Center Community Veloc-
ity Model [
Magistrale et al.
, 2000]. We normalize stress
related fault constitutive model parameters with respect to
the shear modulus, which gives rise to behavior that is
similar to using a homogeneous half-space. We use the
same spontaneous dynamic rupture finite element code as
we have used in previous studies [e.g.,
Aagaard et al.
, 2001,
2004]; in benchmark tests [
Harris and Archuleta
, 2004] it
gives very similar results to several other finite element and
finite difference spontaneous rupture codes. We discretize
the domain using linear tetrahedral finite elements with the
length of element edges set to approximately one tenth of
the wavelength of shear waves with periods of 1.0 s. The
elements at the bottom of the domain have 375 m long
edges whereas the elements at the top of the domain have
265 m long edges (very similar results were obtained for
simulations with elements twice these sizes). This discreti-
zation resolves fault rupture associated with wave propaga-
tion at periods of 1.0 s and longer.
[
11
] Numerous fault constitutive models have been pro-
posed to describe various aspects of slip on faults associated
with earthquake rupture. Some, such as rate and state
friction [e.g.,
Dieterich
, 1979], are derived directly from
laboratory measurements. Others, such as slip weakening
[e.g.,
Ida
, 1972], are used primarily because they are simple
and promote numerical stability. More recently,
Ohnaka
[2003] proposed a slip-based fault constitutive model that
encapsulates both laboratory observations and simple scal-
ing relationships. Alternatively,
Rice
[2006] proposed that
field observations of thin primary slip surfaces and a lack of
widespread melting in fault cores imply thermal processes
(thermal pressurization and flash heating) may control
weakening with the onset of slip in earthquake rupture.
Other fault constitutive models, such as time weakening
[
Andrews
, 2004], are variations of theoretical (in this case,
slip weakening) or laboratory-based models that provide
greater numerical stability for use in spontaneous rupture
simulations. In this study, we use a simple fault constitutive
model that blends the classic slip-weakening friction model
with the restrengthening present in slip-plus-rate-weakening
friction, while including some additional features to ensure
its numerical stability. However, it is important to recognize
that choosing fault constitutive models that promote numer-
ical stability could potentially remove important physics
that occur in real earthquakes.
[
12
] Both slip-weakening and slip-plus-rate-weakening
friction models are widely used fault constitutive models
in spontaneous dynamic rupture modeling [see, e.g.,
Harris
,
2004, and references therein]. In slip-plus-rate-weakening
models friction increases (the fault restrengthens) as slip rate
decreases, whereas in slip-weakening models friction usu-
ally remains at the sliding friction level even after sliding
stops (no restrengthening occurs). In contrast, rate and state
friction models include an increase in friction with the
logarithm of the time once sliding stops. A simple way to
incorporate restrengthening in slip-weakening models is to
impose instantaneous restrengthening wherein the strength
returns to its presliding level immediately upon the termi-
nation of slip [see, e.g.,
Aagaard et al.
, 2001]. However, this
abrupt restrengthening often increases numerical noise in
the simulations through sudden stops and restarts in sliding.
[
13
] In order to generate a continuous spectrum of shear
restrengthening behavior, we develop a common parameter-
ization for slip weakening with and without restrengthening.
Other studies [
Madariaga and Cochard
, 1994;
Shaw
, 1995;
Nielsen et al.
, 2000] have also blended slip-weakening and
slip-plus-rate-weakening friction models using various for-
mulations. We introduce a state variable,
d
, to add restrength-
ening to the linear slip-weakening friction model yielding
s
f
¼
s
fail

s
sliding

1

d
t
ðÞ
D
0

þ
s
sliding
;
_
d
t
ðÞ¼
_
Dt
ðÞ
Ct
ðÞð
1
Þ
Ct
ðÞ¼
(
_
d
r
;
if
_
Dt
ðÞ
V
o
&

Dt
ðÞ
0
0
;
otherwise
_
d
t
ðÞ
_
d
max
;
0

d
t
ðÞ
D
0
;
ð
2
Þ
Table 1.
Description of Parameters
Parameter
Description
s
f
Friction stress
s
fail
Failure stress
s
sliding
Sliding friction stress
s
0
Initial stress
t
xy
Stress tensor component
xy
l
,
m
Lame’s constants
D
(
t
)
Slip at a point on the fault
_
D
(
t
)
Slip rate at a point on the fault
D
0
Slip-weakening parameter in fault
constitutive model
V
0
Restrengthening threshold in fault
constitutive model
C
(
t
)
Restrengthening parameter in fault
constitutive model
d
(
t
)
State variable in fault constitutive model
_
d
(
t
)
Rate of change in state variable
_
d
max
Maximum rate of change in state variable
_
d
r
Rate of restrengthening in fault
constitutive model
k
Rupture propagation parameter,
see equation (10)
L
Length scale related to rupture width,
see equation (10)
W
b
Breakdown work
k
*
Nominal value of
k
a
dyn
Fraction of perturbation in
k
accounted
for by perturbations in dynamic stress drop
a
fail
Fraction of perturbation in breakdown work
accounted for by perturbations in failure s
tress relative to sliding friction stress
Table 2.
Material Properties
a
Depth (km)
Mass Density
(kg/m
3
)
Dilatational Wave
Speed (km/s)
Shear Wave
Speed (km/s)
0
2600
5.70
3.40
11.0
2800
6.30
3.60
22.0
3000
6.90
3.85
38.0
2300
7.80
4.50
40.0
3300
7.80
4.50
a
Control points describing the piecewise linear variation of material
properties with depth.
B04301
AAGAARD AND HEATON: FAULT FRICTION AND SLIP HETEROGENEITY
3of18
B04301
where
s
f
denotes the friction stress,
s
fail
denotes the
failure stress,
s
sliding
denotes the sliding friction stress,
D
denotes slip,
D
0
denotes the slip-weakening parameter, and
_
d
r
,
V
0
, and
_
d
max
are other parameters (see Table 1). Figure 3
shows a typical trajectory of the friction stress in slip and
slip rate space. We choose to specify the friction explicitly
in terms of stresses and independent of the fault normal
stress, because dynamic changes in friction during slip may
be independent of the overburden pressure as a result of
dynamic processes [see, e.g.,
Brune et al.
, 1993;
Melosh
,
1996;
Tworzydlo and Hamzeh
, 1997;
Brodsky and Kanamori
,
2001].
[
14
] Traditional slip weakening behavior with no
restrengthening is produced by setting
_
d
r
= 0. This results
in
C
= 0 for all
t
and
d
(
t
)=
D
(
t
) until
d
(
t
)=
D
0
after which
it is clipped to
D
0
. We also include elements of time-
weakening friction [
Andrews
, 2004] to increase the numer-
ical stability of the solution by limiting the maximum rate of
decrease of the friction stress,
_
d
(
t
)

_
d
max
. This increases the
size of the breakdown zone to several hundred meters for
several meters of slip. This breakdown zone is 2 or more
orders of magnitude larger than those in natural earthquakes
if processes such as flash heating cause very rapid dynamic
shear stress changes. In fact, rapid changes in dynamic
friction necessarily imply small breakdown work. Thus,
while we use a breakdown work comparable to radiated
energy in our numerical model to promote numerical
stability, we recognize that breakdown work could be much
smaller and that the spatial distribution of prestress may
largely control earthquake dynamics.
[
15
] For the restrengthening portion of the fault constitu-
tive model, we employ a simple dependence on a slip rate
threshold,
V
0
, and a rate of restrengthening,
_
d
r
.If
_
d
r
is
chosen to be nonzero (and positive), then when the slip is
decelerating (

D
(
t
)

0) and the slip rate is below some
threshold (
V
0
),
C
> 0 and the state variable decreases toward
zero. As a result, friction increases associated with
restrengthening. If
V
0
= 0, then restrengthening at a point
occurs only after sliding terminates. This corresponds to
traditional slip-weakening behavior during sliding but the
fault restrengthens at a constant rate when sliding stops; this
is somewhat faster than the logarithmic time dependence
found in rate and state friction [e.g.,
Dieterich
, 1979]. On
the other hand, if
V
0
> 0, then restrengthening occurs while
the slip rate is greater than zero producing rate-activated
restrengthening where
_
d
r
controls the rate of restrengthen-
ing. This formulation is slightly more stable numerically than
a direct rate dependence (used in previous studies such as
those by
Aagaard et al.
[2001]
Anderson et al.
[2003], and
Aagaard et al.
[2004]) because the rate of restrengthening is
constant and does not increase as the slip rate decreases.
[
16
] We consider three levels of restrengthening: no
restrengthening (
_
d
r
= 0), slip weakening with restrength-
ening after sliding stops (
V
0
=0,
_
d
r
= 0.325 m/s), and
slip weakening with rate-activated restrengthening (
V
0
=
0.20 m/s,
_
d
r
= 0.325 m/s). For brevity, we will generally
refer to these three cases as without restrengthening,
restrengthening at rest, and rate-activated restrengthening.
Table 3 summarizes the differences in friction model
parameters for these three levels of restrengthening. In all
three cases, we set
_
d
max
=2.0
D
0
/s so that it takes at least
0.5 s for the friction to drop from
s
fail
to
s
sliding
and we can
resolve the breakdown zone. The remaining two principal
parameters of the fault constitutive model (
s
fail
and
D
0
) will
be included in the spatial variation of the parameters as
discussed later.
Figure 3.
A typical trajectory of the friction stress in slip
(D) and slip rate (V) space. The friction decreases in
response to slip weakening and then increases as the slip
rate drops below a threshold. The gray, dashed portion of
the trajectory depicts the restrengthening after sliding stops.
Figure 2.
(top) Material properties and (bottom) fault
constitutive properties as a function of depth. The shear
wave speed (
v
s
), dilatational wave speed (
v
p
), and mass
density (
r
) define the material properties. The sliding stress
(
s
sliding
) increases linearly with depth, and the failure stress
(
s
fail
) varies with the shear modulus.
B04301
AAGAARD AND HEATON: FAULT FRICTION AND SLIP HETEROGENEITY
4of18
B04301
[
17
] We construct the prestress field by superposition of
(1) a nominal background stress field, (2) a uniform shear
strain field, and (3) a random strain field. The nominal
background stress field,
t
xy
¼
6 MPa

1
:
75 MPa
=
km
z
;
ð
3
Þ
t
xx
¼
1
:
0

10

4
4
ml
þ
m
ðÞ
l
þ
2
m
þ
Z
z
0
r
s
ðÞ

r
w
ðÞ
gds
;
ð
4
Þ
t
yy
¼
1
:
0

10

4
2
ml
l
þ
2
m
þ
Z
z
0
r
s
ðÞ

r
w
ðÞ
gds
;
ð
5
Þ
t
zz
¼
Z
z
0
r
s
ðÞ

r
w
ðÞ
gds
;
ð
6
Þ
t
yz
¼
t
xz
¼
0
;
ð
7
Þ
corresponds to linearly increasing shear stress with depth
(z is positive upward), lithostatic overburden pressure (with
hydrostatic pore pressure), and horizontal compression,
where
r
(
z
) denotes the density as a function of
z
,
r
w
denotes
the density of water,
g
denotes the acceleration of gravity,
and
l
and
m
denote Lame’s constants. The shear stress
increases with depth at a rate such that heat production
associated with a 2 cm thick zone and 2 m of frictional
sliding (excluding breakdown work or fracture energy) at a
depth of 15 km will give rise to temperature changes
less than 1200 K (using a heat capacity of per unit mass of
1000 J/(kg K)). This modest temperature increase is
consistent with the lack of substantial melting found in
some exhumed faults [
Chester and Chester
, 1998].
[
18
] The overburden pressure dominates the fault normal
and fault parallel axial stresses (
t
xx
and
t
yy
) except near the
surface. With compressive strains in the fault normal
direction (

xx
=1.0

10

4
), we keep the fault in compres-
sion even under dynamic loading. As a result, the angle
between the fault and the maximum principal stress direc-
tion is a little more than 45 degrees except at depths less
than about a kilometer where the overburden pressure is not
significantly larger than the fault normal compressive stress.
[
19
] On the basis of the characterization of slip hetero-
geneity by
Mai and Beroza
[2002] with von Karman
distributions, we construct the random strain field by
finding the strain field associated with a random slip field
by solving the static elasticity equation in our model. The
power spectra of the random slip follows a von Karman
distribution,
Pk
ðÞ¼
a
y
a
z
1
þ
k
2
ðÞ
H
þ
1
ðÞ
;
ð
8
Þ
where
k
is the radial wave number (
k
2
=
a
y
2
k
y
2
+
a
z
2
k
z
2
),
a
y
=
20 km,
a
z
= 10 km, and
H
= 0.75. We build the slip
distribution by starting with a normally distributed random
field (zero mean with a variance of 1.0) in the wave number
domain. The magnitude of the amplitude in the wave
number domain is set to match the von Karman distribution.
We prevent inclusion of short-length-scale features that the
finite element model cannot resolve by setting the spectral
components associated with wavelengths shorter than
5.0 km to zero. After transforming to the spatial domain, we
taper the random distribution along the edges of the fault as
shown in Figure 1. This approach is very similar to the one
used by
Andrews and Harris
[2005]. The tapering limits the
seismogenic zone to the top half of the domain and prevents
ruptures from propagating to the lateral edges of the domain
(where they would create edge effects not present in natural
earthquakes). The taper follows an exponential variation,
f
taper
¼
1

1
2
e
3
s
=
s
0
if
s
<
0
1
2
e

3
s
=
s
0
if
s
0
;
8
>
<
>
:
ð
9
Þ
where
s
denotes distance from the center of the taper, and
s
0
is the nominal taper width (5.0 km). The final step involves
scaling the distribution by a factor of 150 to create the slip
distribution used in the static simulation to calculate the
random strain field on the fault surface.
[
20
] We want to produce similar slip time histories from
different combinations of initial shear stress and fault
constitutive parameters. Unfortunately, we do not have a
general prescription for generating the initial conditions and
fault constitutive model parameters required to produce a
given spatial-temporal evolution of slip. Consequently,
using trial and error we selected a set of uniform fault
constitutive parameters and the spatially heterogeneous
stress field described earlier to generate a magnitude 7.2
event that ruptures a sizable fraction of our fault (820 km
2
).
We use these parameters and rupture as our reference from
which we construct other simulations with different initial
conditions and fault constitutive parameters that yield
similar rupture behavior.
[
21
]
Day
[1982] and
Madariaga and Olsen
[2000] found
that a nondimensional parameter,
k
¼
s
0

s
sliding

2
L
W
b
m
ð
10
Þ
(
s
0
is the initial shear stress,
s
sliding
is the sliding friction,
W
b
is the breakdown work,
m
is the shear modulus, and
L
is
the width of the slip pulse), correlates with the extent and
speed of rupture propagation. We follow
Tinti et al.
[2005]
in using the term breakdown work instead of fracture
energy, because we consider fault constitutive behavior at a
macroscopic scale and do not distinguish among the various
Table 3.
Friction Parameters for the Three Levels of Shear
Restrengthening
a
Scenario
_
d
r
(m/s)
V
0
(m/s)
noR
0
0
restR
0.325
0
rateR
0.325
0.20
a
The friction models noR, restR, and rateR refer to slip weakening
without restrengthening, slip weakening with restrengthening at rest, and
slip weakening with rate-activated restrengthening, respectively.
B04301
AAGAARD AND HEATON: FAULT FRICTION AND SLIP HETEROGENEITY
5of18
B04301
mechanisms dissipating energy at the rupture front, e.g.,
creation of new surfaces (fracture energy), frictional sliding
on secondary slip surfaces, and inelastic deformation in the
surrounding volume. This nondimensional parameter is
proportional to the inverse of the critical crack length
required for rupture propagation [
Andrews
, 1976b;
Day
,
1982] and roughly approximates the ratio of the available
strain energy (
W

(
s
0

s
sliding
)
2
L
/
m
) to the breakdown
work (
W
b
)[
Madariaga and Olsen
, 2000]. If we replace the
sliding friction with the final shear stress in the expression
for
k
this would be less of an approximation but we would
not be able to compute
k
prior to rupture. Although
k
is
formulated in the context of a homogeneous half-space with
uniform initial shear stresses and fault constitutive model
parameters, it provides a convenient parameter related to
macroscopic rupture behavior (rupture extent and rupture
speed) that we can keep constant across our suite of
simulations. Other parameters, such as dynamic stress drop
and relative strength
S
=(
s
fail

s
0
)/(
s
0

s
sliding
)
[
Andrews
, 1976a]), display a greater dependence on the
details of the fault constitutive model for a given spatial-
temporal evolution of slip than
k
.
[
22
] We construct the common baseline for our suite of
spontaneous rupture simulations by computing
k
on the
fault surface for our reference magnitude 7.2 event. We also
compute a nominal, uniform
k
*, such that
k
/
k
* describes
the perturbation in the initial conditions and fault constitu-
tive parameters from some nominal, uniform state. If we
arbitrarily select
L
= 1 km (because
L
loses its original ties
to rupture width in our heterogeneous setting), then
k
= 0.19
corresponds to our nominal, uniform value; we denote the
associated nominal values for the initial shear stress, sliding
friction, and breakdown work by
s
*
0
,
s
*
sliding
,and
W
*
b
,
where
k
*=(
s
*
0

s
*
sliding
)
2
/(
W
*
b
m
) as expected. Table 4
gives the nominal values of the parameters associated with
k
* in our simulations. This nominal value,
k
*, is the
‘‘average’’ over the entire fault surface; the ruptures will
preferentially propagate over the area in which
k
is locally
high, so the average over the area where slip occurs will be
considerably higher than 0.19.
[
23
] We construct our suite of spontaneous rupture sim-
ulations with variations in the spatial distributions of the
initial shear stress and fault constitutive parameters by
associating these variations in the parameters with the
spatial variation in
k
(see Figures 4–6 top) relative to
k
*.
This provides a systematic way of considering a continuous
range of trade-offs among the dynamic stress drop (initial
shear stress) and fault constitutive parameters (e.g., break-
down work, slip-weakening parameter, and failure stress).
[
24
] The random strain field used to construct the spatial
variation in
k
contains nonzero tractions in both the along-
strike and updip directions. However,
k
is a scalar, so it
does not contain information about the perturbations in the
orientation of the initial stress field; consequently, we
simply apply the shear stresses in the horizontal direction,
which is consistent with the strike-slip faulting. We allow
rake rotations so the slip is not constrained to be in the
horizontal direction.
[
25
] Examining the expression for
k
, one recognizes that
mapping perturbations in
k
results in examining trade-offs
between perturbations in dynamic stress drop (
s
0

s
sliding
)
and breakdown work (
W
b
). In our parameterization break-
down work is characterized by the failure stress (
s
fail

s
sliding
) and the slip-weakening parameter (
D
0
). Appendix A
explains how we relate perturbations in
k
to perturbations in
the initial shear stress, failure stress, and slip-weakening
parameter using nondimensional parameters
a
dyn
and
a
fail
.
a
dyn
specifies the fraction of the perturbation in
k
that arises
from perturbations in the dynamic stress drop, and
a
fail
specifies the fraction of the perturbation in the breakdown
work that arises from perturbations in the failure stress. This
yields
s
0

s
sliding
m
¼
s
0
*

s
sliding
*
m

k
k
*
a
dyn
2
ð
11
Þ
for the initial stress field relative to the sliding stress
(dynamic stress drop),
s
fail

s
sliding
m
¼
s
fail
*

s
sliding
*
m

k
k
*
a
fail
a
dyn

1
ðÞ
ð
12
Þ
for the failure stress relative to the sliding friction stress, and
D
0
¼
D
0
*
k
k
*
1

a
fail
ðÞ
a
dyn

1
ðÞ
ð
13
Þ
for the slip-weakening parameter, where
a
dyn
and
a
fail
can
be between 0 and 1, inclusively. When
a
dyn
= 1, spatial
variations in
k
arise exclusively from spatial variations in
dynamic stress drop; alternatively, when
a
dyn
= 0, spatial
variations in
k
arise exclusively from spatial variations in
breakdown work. Similarly, when
a
fail
=1,spatial
variations in breakdown work arise exclusively from spatial
variations in the failure stress relative to the sliding friction
stress, whereas when
a
fail
= 0, spatial variations in
breakdown work arise exclusively from spatial variations
in the slip-weakening parameter,
D
o
.
[
26
] In this study we focus on the spatial variability of the
initial and final shear stress relative to the sliding friction
stress, rather than the spatial variability in the absolute level
of the stress. For our sets of parameters, including spatial
variations in the sliding friction, while keeping the initial
and failure stress relative to the sliding friction the same,
does not change the spatial and temporal evolution of slip.
Furthermore, we will use these relative stresses to judge the
preearthquake and postearthquake consistency in the spatial
heterogeneity of the stress field. Consequently, we do not
consider spatial variations in the sliding friction: we map
spatial variations in the dynamic stress drop into spatial
variations in the initial shear stress, and spatial variations in
the failure stress relative to the sliding friction directly into
the failure stress, keeping the sliding friction stress equal to
Table 4.
Parameters Associated With Nominal Value of
k
a
L
*
1.00 km
W
*
b
/
m
1.62

10

4
D
*
0
0.65 m
(
s
*
fail

s
*
sliding
)/
m
5.00

10

4
(
s
*
0

s
*
sliding
)/
m
1.76

10

4
a
With
k
* = 0.19. Note that the nominal sliding stress
s
*
sliding
is equal to
the nominal background stress
t
xy
(equation (3)).
B04301
AAGAARD AND HEATON: FAULT FRICTION AND SLIP HETEROGENEITY
6of18
B04301
Figure 4.
Comparison of initial conditions and rupture characteristics for three illustrative scenarios for
slip-weakening friction without restrengthening. All scenarios use the same initial distribution of
k
(top
middle). The spatial variation in
k
is associated with either variations in dynamic stress drop via the
prestress (
a
dyn
= 1; left), variations in dynamic stress drop via the prestress accompanied by variations in
breakdown work via
D
o
(
a
dyn
= 0.75,
a
fail
= 0; middle), or variations in dynamic stress drop via the
prestress accompanied by variations in breakdown work via both
D
o
and failure stress (
a
dyn
= 0.75,
a
fail
=
0.25; right).
B04301
AAGAARD AND HEATON: FAULT FRICTION AND SLIP HETEROGENEITY
7of18
B04301
Figure 5.
Comparison of initial conditions and rupture characteristics for three illustrative scenarios for
slip-weakening friction with restrengthening at rest. All scenarios use the same initial distribution of
k
(top middle). The spatial variation in
k
is associated with either variations in dynamic stress drop via the
prestress (
a
dyn
= 1; left), variations in dynamic stress drop via the prestress accompanied by variations in
breakdown work via
D
o
(
a
dyn
= 0.75,
a
fail
= 0; middle), or variations in dynamic stress drop via the
prestress accompanied by variations in breakdown work via both
D
o
and failure stress (
a
dyn
= 0.75,
a
fail
=
0.25; right).
B04301
AAGAARD AND HEATON: FAULT FRICTION AND SLIP HETEROGENEITY
8of18
B04301
Figure 6.
Comparison of initial conditions and rupture characteristics for three illustrative scenarios for
slip-weakening friction with rate-activated restrengthening. All scenarios use the same initial distribution
of
k
(top middle). The spatial variation in
k
is associated with either variations in dynamic stress drop via
the prestress (
a
dyn
= 1; left), variations in dynamic stress drop via the prestress accompanied by variations
in breakdown work via
D
o
(
a
dyn
= 0.75,
a
fail
= 0; middle), or variations in dynamic stress drop via the
prestress accompanied by variations in breakdown work via both
D
o
and failure stress (
a
dyn
= 0.75,
a
fail
=
0.25; right).
B04301
AAGAARD AND HEATON: FAULT FRICTION AND SLIP HETEROGENEITY
9of18
B04301
the nominal background stress. On the other hand, in real
earthquakes the sliding friction likely does vary spatially
and its absolute level and spatial variability become impor-
tant in controlling slip and stress heterogeneity as inelastic
deformation, heating, and other processes occur.
[
27
] The two trade-off parameters,
a
dyn
and
a
fail
, permit
selecting from a continuous spectrum of initial conditions
and fault constitutive parameters that match the selected
spatial variation in
k
. We consider four values for
a
dyn
and
three values for
a
fail
covering a considerable range of the
parameter space; Table 5 gives the values of the trade-off
parameters for the 10 cases in this study. We consider these
10 cases for each of the three levels of shear restrengthen-
ing, yielding a total of 30 scenarios.
3. Results
[
28
] In each scenario the ruptures nucleate at a depth of
8.3 km at a location 43.4 km along strike from the left (

y
)
edge of the fault as a result of raising the shear stress to 2%
above the failure stress over a circular region with a 2 km
radius. This region is centered on the location where
k
is at
a maximum. Table 5 summarizes the basic characteristics of
fault rupture in each of the 30 scenarios, while Figures 4–6
illustrate examples of the rupture behavior and conditions
before and after rupture for the three levels of shear
restrengthening. We succeed in producing 16 ruptures with
similar rupture behavior: bilateral rupture with moment
magnitudes in the range of 6.8 to 7.2 and similar spatial-
temporal evolutions of slip.
[
29
] Figure 7 shows rupture propagation via snapshots
of the slip rate 6.0 s after the ruptures begin, and Figure 8
demonstrates that the velocity time histories from these
scenarios are quite similar, and in some cases nearly identical.
It also shows how delayed rupture to the left (

y
direction) in
a few of the scenarios affects the ground motions, indicating
that such deviations in rupture behavior could be identified in
a kinematic source inversion. The scenarios labeled ‘‘dis-
carded’’ in Table 5 produce ruptures with significantly
different behavior suggesting that other factors besides the
value of
k
influence the extent of rupture. For example, the
ratio of the failure stress relative to the dynamic stress drop
and the initial stress, as parameterized by the relative
strength,
S
=(
s
fail

s
0
)/(
s
0

s
sliding
)[
Andrews
, 1976a],
is another nondimensional parameter that influences the
ability of ruptures to propagate, particularly in discrete
models that blunt the sharpness of the stress peak at the
leading edge of the rupture.
Table 5.
Parameters and Rupture Characteristics for Each Scenario
a
a
dyn
a
fail
Average Slip (m)
M
w
SD in Strength Excess,
b

10

4
Remarks
Before
After
Friction Model noR
1.00
2.81
7.19
1.2
0.2
Inconsistent
0.75
0.00
1.99
7.04
1.0
0.2
Inconsistent
0.75
0.25
1.91
7.02
1.0
0.6
Marginally self-consistent
0.75
0.50
1.86
7.01
1.1
0.7
Marginally self-consistent
0.50
0.00
1.61
6.93
0.7
0.2
Inconsistent
0.50
0.25
1.51
6.89
1.0
0.7
Marginally self-consistent
0.50
0.50
1.44
6.86
1.2
0.8
Marginally self-consistent
0.25
0.00
1.34
6.85
0.5
0.2
Inconsistent
0.25
0.25
1.13
6.64
Small event (discarded)
0.25
0.50
0.63
6.23
Small event (discarded)
Friction Model restR
1.00
2.62
7.16
1.2
0.4
Inconsistent
0.75
0.00
1.81
7.01
0.9
0.5
Inconsistent
0.75
0.25
1.77
7.00
1.0
0.6
Marginally self-consistent
0.75
0.50
1.72
6.98
1.1
0.7
Marginally self-consistent
0.50
0.00
1.31
6.84
0.7
0.6
Self-consistent
0.50
0.25
1.22
6.66
Small event (discarded)
0.50
0.50
1.15
6.63
Small event (discarded)
0.25
0.00
1.04
6.59
Small event (discarded)
0.25
0.25
0.93
6.54
Small event (discarded)
0.25
0.50
0.60
6.20
Small event (discarded)
Friction Model rateR
1.00
2.35
7.12
1.2
0.6
Inconsistent
0.75
0.00
1.30
6.84
0.9
0.8
Self-consistent
0.75
0.25
1.29
6.83
1.0
0.7
Marginally self-consistent
0.75
0.50
1.23
6.68
Small event (discarded)
0.50
0.00
0.99
6.58
Small event (discarded)
0.50
0.25
0.94
6.56
Small event (discarded)
0.50
0.50
0.71
6.29
Small event (discarded)
0.25
0.00
0.71
6.25
Small event (discarded)
0.25
0.25
0.66
6.21
Small event (discarded)
0.25
0.50
0.52
6.10
Small event (discarded)
a
The friction models noR, restR, and rateR refer to slip-weakening without restrengthening, slip-weakening with restrengthening at rest, and slip-
weakening with rate-activated restrengthening, respectively.
b
The values for the standard deviation in strength excess have been normalized by the shear modulus.
B04301
AAGAARD AND HEATON: FAULT FRICTION AND SLIP HETEROGENEITY
10 of 18
B04301
[
30
] Rupture propagation beginning at the hypocenter and
proceeding to the left (

y
direction) is more difficult than
that to the right (+
y
direction) due to a large local minimum
in
k
just to the left of the hypocenter. As more of the spatial
variation in
k
is mapped into spatial variations in break-
down work and failure stress, small regions with locally
large values of strength excess prove to be significant
barriers for propagation in the discrete models.
k
is inversely
proportional to the breakdown work (and therefore the
failure stress) but proportional to the square of the differ-
ence between initial shear stress and sliding friction. This
means that variations in the dynamic stress drop (initial
shear stress) result in relatively larger variations in
k
than do
variations in failure stress. Likewise, variations in
k
arise
from smaller variations in the dynamic stress drop (initial
shear stress) compared with the failure stress (compare the
red surfaces in Figures 4 (middle), 5 (middle), and 6
(middle)). As a result, mapping a significant portion of
the spatial variation in
k
into the failure stress (
a
dyn

0.5)
leads to dramatic variations in the failure stress preventing
the ruptures from expanding beyond the nucleation region.
This is why scenarios with
a
dyn

0.5 produce earthquakes
(in all but a couple cases) with average slips and moment
magnitudes that are smaller than the target values. Addi-
tionally, the presence of shear restrengthening, which
increases the friction as the slip rate becomes low or when
sliding stops, tends to extinguish the rupture if it encounters
any sort of barrier (locally small
k
). This explains why the
events tend to become smaller as more restrengthening is
added to the fault constitutive model; rapid restrengthening
promotes slip heterogeneity, which also increases the pos-
sibility of rupture termination [
Liu-Zeng et al.
, 2005].
[
31
] The presence of restrengthening permits the dynamic
stress drop to be significantly larger than the static stress
drop, but the difference between the dynamic and static
stress drop is very small over most of the rupture area. This
indicates that just a small rise in the dynamic sliding friction
during slip can have a significant affect on the healing
phases (i.e., pulse width) and whether or not the rupture
propagates. However, for the pulse widths permitted with
the spatial resolution of our models, it does not have a
strong affect on the macroscopic rupture behavior and the
conditions left behind by the rupture, aside from adding in a
small amount of short-length-scale numerical noise. On the
other hand, the spatial variations in the breakdown work
(manifested in variations in the failure stress and slip-
weakening parameter) create barriers that impede slip,
resulting in final shear stress that differs significantly at
many locations from the sliding stress. These barriers can
lead to slight variations in rupture speed and slip rate as is
evident in the more complex velocity time histories at sites
A and B in Figure 8 for the cases with
a
dyn

0.75,
nevertheless the ground motions remain quite similar.
[
32
] We now focus our attention on the initial and final
shear stress relative to the sliding friction and evaluate the
consistency in the amount of heterogeneity. Our objective is
to identify the sets of parameters conducive to generating
‘‘self-consistent’’ ruptures; as described earlier, we perceive
self-consistent ruptures as those that begin with a heteroge-
neous shear stress field and leave behind a heterogeneous
shear stress field. We quantify the degree of self-consistency
by comparing the standard deviation in the normalized
strength excess before and after rupture. At each point
where slip occurs, we compute the strength excess (the
Figure 7.
Snapshots of slip rate 6.0 s after rupture initiation for each of the nine scenarios in
Figures 4–6, left to right corresponding to the same three different realizations of initial conditions and
parameters shown in Figures 4–6 and the top to bottom corresponding to the three levels of
restrengthening. Whereas the leading edges of the ruptures are very similar, the widths of the slip pulses
differ significantly.
B04301
AAGAARD AND HEATON: FAULT FRICTION AND SLIP HETEROGENEITY
11 of 18
B04301
difference between the failure stress and the shear stress)
and divide by the shear modulus to find the normalized
strength excess. Before rupture, the distribution of normal-
ized strength excess roughly matches a normal distribution,
so it is straightforward to compute the standard deviation for
each of the scenarios. The normalized strength excess after
rupture is bimodal: a peak near zero corresponds to the edge
of the rupture area, whereas a second peak corresponds to
the area where the stress decreases as a result of the rupture.
We focus on this second peak because the standard devia-
tion with respect to this peak gives an indication of the
heterogeneity in the stress field where rupture occurred. If
the standard deviation decreases dramatically, then it indi-
cates that the rupture removed heterogeneity in the shear
stress field; we label these cases ‘‘inconsistent’’. On the
other hand, if the standard deviation is about the same
before and after rupture, then it indicates that the rupture
maintained heterogeneity in the shear stress field; we label
these cases ‘‘self-consistent’’. We label intermediate cases
‘‘marginally self-consistent.’’
[
33
] Figures 4–6 show the strength excess before and
after rupture for nine of the scenarios, and Figure 9 displays
the histograms of normalized strength excess over the
rupture area for each of these scenarios. Table 5 lists the
normalized strength excess of the rupture area for all 16 of
the ruptures near the target size. The standard deviation of
the normalized strength excess differs slightly from one
scenario to another because the rupture areas differ.
Figures 4 (left), 5 (left, 6 (left), and 9 (left) illustrate cases
in which the initial shear stress is heterogeneous but the
final shear stress field is homogeneous over the area where
rupture occurred. As a result, the standard deviation in the
normalized strength excess after rupture is very small, and
we classify these three simulations as ‘‘inconsistent.’’
Figures 4 (middle and right), 5 (middle and right), 6 (middle
and right), and 9 (middle and right) demonstrate that as we
distribute more of the spatial variation in
k
into the
breakdown work, the standard deviations in the normalized
strength excess do not decrease as much, so we begin to
retain heterogeneity in the stress field. Additionally, as we
add restrengthening to the fault constitutive model and
include rate-activated restrengthening, the ruptures tend to
become more self-consistent. Furthermore, other realiza-
tions of spatially heterogeneous
k
constructed in the same
manner produce similar results. Including restrengthening
also results in more discarded simulations for the cases in
which we map at least 50% of the perturbations in
k
into the
breakdown work (giving rise to very large spatial variations
in the failure stress and slip-weakening parameter) as
discussed earlier.
Figure 8.
Velocity time histories of the fault-perpendicular
particle motion at two sites on the ground surface (see
Figure 1) for the 16 scenarios with
M
w
6.8 ruptures (see
Table 5). The labels indicate the friction model and the
values for
a
dyn
and
a
fail
. The models with
a
dyn
0.75
produce very similar rupture behavior and waveforms. The
other scenarios have similar final slip distributions but take
longer to break through the rupture barrier to the left (

y
direction), resulting in variations in waveform shape and
timing.
B04301
AAGAARD AND HEATON: FAULT FRICTION AND SLIP HETEROGENEITY
12 of 18
B04301
[
34
] In all of the simulations slip is concentrated near the
leading edge of the rupture front as shown in Figure 7.
However, the total width of the slip pulses differs depending
on the sources of heterogeneity (comparing columns) and
the amount of restrengthening in the fault constitutive
model (comparing rows). This means that the risetimes
associated with the majority of the slip at a point are about
the same, but the risetimes associated with the final slip
values differ. Kinematic source inversions would not be able
to detect these subtle differences, because the resulting
ground motions are so similar. Incorporating spatial hetero-
geneity in both the dynamic stress drop (initial shear stress)
and breakdown work trims the trailing edge of the slip pulse
as does adding restrengthening to the fault constitutive
model. One might expect the slip-weakening with rate-
activated restrengthening model to consistently produce
the narrowest slip pulses. Instead, the pulse width exhibits
a dependence on both the spatial variation in the parameters
and the presence of restrengthening in the fault constitutive
model. Combining spatial variations in dynamic stress drop
Figure 9.
Histograms of normalized strength excess over the rupture area before and after rupture for
each of the nine scenarios in Figures 4–6, with left to right corresponding to the same three different
realizations of initial conditions and parameters shown in Figures 4–6 and top to bottom corresponding
to the three levels of restrengthening. Including spatial heterogeneity in both the initial shear stress and
breakdown work (especially in the failure stress) results in only slight decreases in the standard deviation,
indicating the ruptures retain most of the heterogeneity in the stress field as desired.
B04301
AAGAARD AND HEATON: FAULT FRICTION AND SLIP HETEROGENEITY
13 of 18
B04301
and breakdown work aid the generation of healing fronts,
contributing to narrow slip pulses.
4. Implications
[
35
] There are two important contexts in which to exam-
ine the results of these simulations. First, the simulations
suggest some useful constraints on fault constitutive models
for use in numerical simulations of earthquake ruptures that
focus on reproducing realistic distributions of both slip and
stress. Second, we consider how to extrapolate the numer-
ical results, which span a limited range of scales, to the
dynamics of natural earthquakes, which span a much
broader range of scales.
4.1. Creating Realistic Numerical Simulations
[
36
] Our results imply that generating earthquake ruptures
in numerical models consistent with the basic physics
driving heterogeneous slip requires spatially heterogeneous
fields for both the dynamic stress drop and breakdown
work. Models with only heterogeneous dynamic stress drop
can produce rupture behavior similar to those that include
spatial heterogeneity in both, but the initial and final stress
fields are not self-consistent. For example, the shear stresses
in the scenarios with heterogeneous initial shear stresses and
homogeneous breakdown work became homogeneous in
the area where rupture occurred. Such models might be
acceptable for use in studying ground motions, but they
neglect complex spatial variations in the initial conditions
and parameters that appear to play important roles in the
physics of earthquake ruptures. As a result, we consider
such models to be inconsistent with the conditions driving
heterogeneous slip in earthquakes.
[
37
] In this study we focused on the spatial variability of
the initial and final shear stresses relative to the sliding
friction stress and considered spatial variations in dynamic
stress drop arising from spatial variations in the initial shear
stress alone. In the broader context of real earthquakes, we
need to consider how such spatial variations in the dynamic
stress drop might develop from spatial variations arising in
either the shear stresses or the sliding friction. Seismologic
observations support characterization of the stress field as
heterogeneous both before and after rupture. The rapid
onset of aftershocks, many of which appear to coincide
with the main shock rupture surfaces [see, e.g.,
Bakun et
al.
, 2005;
Hauksson and Shearer
, 2005;
Langbein et al.
,
2005], suggests that the strength excess following an
earthquake rupture is quite heterogeneous. Associating this
heterogeneous strength excess with heterogeneity in the
stress field would explain the presence of left-lateral after-
shocks following the 1989 Loma Prieta earthquake, which
had oblique right-lateral slip [
Oppenheimer
, 1990]. Stress
tensor orientations inferred from focal mechanism inversions
suggest that, in general, the stress field is spatially hetero-
geneous [
Hardebeck and Hauksson
, 2001a, 2001b]. Further-
more,
Smith
[2006, available at http://resolver.caltech.edu/
CaltechETD:etd-05252006-191203] used a simple 3-D
model of a heterogeneous stress field to explain the rela-
tionship among tectonic loading, stress perturbations due to
earthquake rupture, and temporal variations in focal mech-
anisms. Stronger evidence for stress heterogeneity comes
from measurements of crustal stress from the Cajon Pass
and KTB boreholes, which support characterization of the
crustal stress field as a fractal distribution [
Marsan and
Bean
, 2003].
[
38
] We do not have similar evidence for spatial hetero-
geneity in the sliding friction or breakdown work, because
in situ measurements of spatial variations of these quantities
are, unfortunately, not currently practical.
Madariaga
[1983] suggested that heterogeneous breakdown work could
explain high-frequency radiation. While the analysis was
limited to a semi-infinite crack, we expect heterogeneous
breakdown work has a similar affect on narrow slip pulses.
Several studies [
Guatteri et al.
, 2003;
Ide
, 2003;
Tinti et al.
,
2005] have made indirect estimates of the spatial variations
in breakdown work using the slip-weakening model and
seismologic data. Additional indirect evidence for the pos-
sible existence of spatial variations in sliding friction and
breakdown work comes from studies of complex fault
geometry, such as irregular nonplanar surfaces, bends, and
junctions [
Segall and Pollard
, 1980;
Andrews
, 1989;
Saucier
et al.
, 1992;
Harris and Day
, 1993;
Bouchon
, 1997;
Nielsen
and Knopoff
, 1998]. It is relatively straightforward to include
nonplanar fault geometry at length scales much larger than
the discretization size in numerical simulations, but it is not
clear how to implement geometrically complex fault surfaces
at small scales or the physics associated with slip on such
surfaces [
Harris
, 2004]. Recent attempts to develop equiva-
lent planar representations of kinks or bends with fault
surfaces have had marginal success [
Kase and Day
, 2006;
Madariaga et al.
, 2006].
[
39
] In our simulations we find that restrengthening may
not be an essential ingredient for producing self-consistent
behavior, but it does expand the range of conditions under
which ruptures exhibit self-consistent behavior. The most
important ingredient for developing self-consistent ruptures
appears to be narrow slip pulses. Including spatial hetero-
geneity in the breakdown work in addition to the dynamic
stress drop produces narrow slip pulses, independent of the
presence or absence of restrengthening. This is not unex-
pected given that several previous studies have shown that
heterogeneous stress fields and rate-independent friction are
compatible with the pulse-like behavior imaged in kinematic
source inversions [
Das and Aki
, 1977;
Papagerogiou and
Aki
, 1983;
Beroza and Mikumo
, 1996;
Ide and Takeo
, 1997;
Day et al.
, 1998;
Fukuyama and Madariaga
, 1998;
Oglesby
and Day
, 2002].
[
40
] In this study we use a rather simple test for gauging
whether the ruptures are self-consistent. A more rigorous
test could involve modeling multiple earthquake cycles that
include strain accumulation, strain release through the
propagation of ruptures, and postseismic relaxation. In such
simulations one can examine long-term trends and quantify
any significant changes in the spectral character of the
distributions of slip over time. This approach is attractive
because several studies [e.g.,
Somerville et al.
, 1999;
Mai
and Beroza
, 2002;
Lavalle
́e and Archuleta
, 2003;
Liu-Zeng
et al.
, 2005] have focused on characterizing the spectral
character of slip distributions. Additionally, aftershocks,
events on adjacent portions of the fault, and nonuniform
loading, which might all influence the evolution of the
stress field, could also be included.
[
41
] Some work with multicycle earthquake simulations
has been done for simplistic 2-D [e.g.,
Shaw
, 1997;
Lapusta
B04301
AAGAARD AND HEATON: FAULT FRICTION AND SLIP HETEROGENEITY
14 of 18
B04301