Close
About
FAQ
Home
Collections
Login
USC Login
Register
0
Selected
Invert selection
Deselect all
Deselect all
Click here to refresh results
Click here to refresh results
USC
/
Digital Library
/
University of Southern California Dissertations and Theses
/
Dynamic rupture processes and seismic radiation in models of earthquake faults separating similar and dissimilar solids
(USC Thesis Other)
Dynamic rupture processes and seismic radiation in models of earthquake faults separating similar and dissimilar solids
PDF
Download
Share
Open document
Flip pages
Contact Us
Contact Us
Copy asset link
Request this asset
Transcript (if available)
Content
DYNAMIC RUPTURE PROCESSES AND SEISMIC RADIATION IN MODELS OF
EARTHQUAKE FAULTS SEPARATING SIMILAR AND DISSIMIAR SOLIDS
by
Zheqiang Shi
A Dissertation Presented to the
FACULTY OF THE GRADUATE SCHOOL
UNIVERSITY OF SOUTHERN CALIFORNIA
In Partial Fulfillment of the
Requirements for the Degree
DOCTOR OF PHILOSOPHY
(GEOLOGICAL SCIENCES)
August 2009
Copyright 2009 Zheqiang Shi
ii
Acknowledgements
First and foremost, I would like to express sincere gratitude to my advisor, Professor
Yehuda Ben-Zion, for his mentoring, support and guidance during every phase of my
Ph.D. career. I thank Yehuda for sharing with me his vast knowledge, deep insight and
critical thinking skills in science.
I am deeply grateful to Professor Alan Needleman, without whose immense help this
thesis would only be half complete. I thank Alan for so much I learned from him during
our collaboration. His enthusiasm and always-positive attitude in doing scientific re-
search are tremendously inspiring.
I would like to thank the members of my thesis committee, who are also on my quali-
fication exam committee, Professors Charles Sammis, Aiichiro Nakano, Thomas Jordan
and Thorsten Becker for their support and advice.I am also very thankful to the depart-
ment staff members, especially Cindy Waite and John McRaney, who helped me so much
in various ways. Special thanks goes to all of my friends and colleagues. Po Chen, now Professor
Chen, provided me much needed support in my first couple of years at USC and has been
a great friend ever since. I learned freestyle swimming, a lot of American culture and
many other things from Jeremy Zechar, who is always a great company to have around.
Youlin Chen and Yufang Rong are like my elder brother and sister who took care of me
iii
every time I went to Boston. Professor Zhigang Peng provided comfortable accommoda-
tion during my two weeklong visits at Georgia Tech and shared with me his expertise on
seismic data processing and valuable personal experience in academia. Finally, I am forever indebted to my family, especially my mother and my sister for
their patience, understanding and endless unconditional support. I also want to thank my
dear aunts and uncles for their enormous caring for so many years.
iv
Table of Contents
Acknowledgements ii
List of Tables vi
List of Figures vii
Abstract xv
Chapter 1 Introduction 1
Chapter 2 Dynamic Rupture Along a Bimaterial Interface Governed
by Slip-Weakening Friction
6
2.1 Introduction 6
2.2 Analysis 10
2.2.1 Model setup 10
2.2.2 A homogeneous case 14
2.2.3 Corresponding cases with material contrast 15
2.2.4 Variable differences between static friction and initial
shear stress
20
2.2.5 Variable differences between the static and dynamic coef-
ficients of friction
23
2.2.6 Supershear ruptures with a larger V
force
29
2.3 Discussion 31
Chapter 3 Properties of Dynamic Rupture and Energy Partition in a
Solid with a Frictional Interface
36
3.1 Introduction 36
3.2 Model setup 39
3.2.1 Finite element formulation 39
3.2.2 Model configurations 40
3.2.3 Constitutive relations 43
3.2.4 Nucleation procedure 48
3.2.5 Energy calculations 48
3.3 Simulation results 50
3.4 Discussion 63
Chapter 4 Slip modes and partitioning of energy during dynamic
frictional sliding between identical elastic-viscoplastic sol-
ids
68
4.1 Introduction 68
4.2 Problem formulation 71
4.2.1 Boundary value problem formulation 71
v
4.2.2 Bulk elastic-viscoplastic constitutive relation 73
4.2.3 Interface constitutive relations 75
4.2.4 Energy partitioning 77
4.3 Numerical results 78
4.3.1 Case I: a crack-like mode 79
4.3.2 Case II: a pulse-like mode 84
4.3.3 Case III: a train-of-pulses mode 89
4.3.4 Parameter studies 93
4.3.4.1 Effect of interface elasticity 93
4.3.4.2 Effect of impact velocity 95
4.4 Discussion 98
4.5 Appendix: a rate tangent method for updating the interface shear
traction
102
Chapter 5 Seismic Radiation From Tensile and Shear point Disloca-
tions Between Similar and Dissimilar Solids
104
5.1 Introduction 104
5.2 Analysis 107
5.2.1 Radiation patterns from point sources in a homogeneous
isotropic solid
107
5.2.2 Synthetic seismograms for dislocation sources in a homo-
geneous solid
109
110
5.2.3 Synthetic seismograms for dislocation sources between
dissimilar solids
117
5.3 Discussion 126
5.4 Appendix: derivation of radiation patterns from shear and tensile
point-dislocation sources
128
Bibliography 133
vi
List of Tables
2.1 Parameter values that are common to all calculations. 13
2.2 Rupture velocities for cases with different parameters. 18
vii
List of Figures
2.1 (a) A Model for 2D finite-difference simulations of mode-II dy-
namic rupture along a bimaterial interface. The initial normal and
shear stress components at the remote boundaries are −σ
∞
and τ
∞
.
Events are nucleated by in two limited x−t domains that expand
symmetrically about the origin. (b) A linear slip-weakening fric-
tion law that governs the propagation of spontaneous rupture out-
side the source region.
11
2.2 (a) Profiles of slip (top panel) and slip velocity (bottom panel)
along the fault at different times for a homogeneous case. (b) The
corresponding spatio-temporal evolution of slip velocity. The re-
sults are symmetric in the positive and negative directions.
14
2.3 (a) Profiles of slip (top panel) and slip velocity (bottom panel)
along the fault at different times for a case with 4% material con-
trast. (b) The corresponding spatio-temporal evolution of the slip
velocity.
16
2.4 Profiles of slip and slip velocity along the fault at different times
for a case with 5% material contrast across the interface. The rup-
ture is arrested quickly in the negative direction and becomes self-
sustained in the positive direction.
17
2.5 Profiles of slip along the fault at different times for cases with
10% (top) and 20% (bottom) material contrast across the interface.
The ruptures propagate in the positive direction in a self-sustaining
manner with increasing slip.
18
2.6 Maximum propagation distances of ruptures in the opposite direc-
tions vs. the degree of material contrast. Cases with self-sustained
rupture in the positive direction are marked with “Divergence”.
The inset gives the ratios of the maximum propagation distance of
rupture in the positive direction divided by the corresponding dis-
tance in the negative direction.
19
2.7 Scalar potencies of ruptures in the opposite directions versus de-
grees of material contrast for several different values of initial
shear stress. The inset shows the ratios of rupture potency in the
positive direction divided by the corresponding value in the nega-
tive directions.
21
viii
2.8 A phase diagram illustrating the effects of µ
s
σ
∞
− τ
∞
and degree of
material contrast on symmetry properties of ruptures along a bima-
terial interface. The limit values of 0 (deep blue) and 2 (dark red)
represent, respectively, symmetric bilateral and fully asymmetric
unilateral ruptures. Detailed results associated with parameter val-
ues marked by the white dots on the black horizontal line are
shown in figures with numbers indicated in the circles. The label
‘5t’ stands for the top panel of Fig. 2.5.
22
2.9 A phase diagram illustrating the effects of µ
s
– µ
d
and degree of
material contrast on symmetry properties of ruptures along a bima-
terial interface. The other symbols and notations are the same as in
Fig. 2.8.
24
2.10 Profiles of slip velocity along the fault at different times for 20%
material contrast across the interface and different values of µ
s
–
µ
d
. The high-frequency oscillations near µ
s
– µ
d
= 0.4 are associ-
ated with a transition from unilateral wrinkle-like rupture pulses
(a-b) to bilateral crack-like ruptures (e-f).
26
2.11 Spatio-temporal evolution of slip velocity for a case with µ
s
– µ
d
=
0.7. Supershear crack fronts propagate in both directions with ve-
locities close to the P-wave speed. The crack front in the positive
direction is followed by a slip pulse with a velocity close to c
GR
.
28
2.12 A summary diagram extended from Fig. 2.9. The dashed line indi-
cates a tentative transition that separates dynamic regimes associ-
ated with predominantly bilateral and predominantly unilateral
ruptures.
29
2.13 Spatio-temporal evolution of slip velocity for a case with 20%
contrast, initial shear stress τ
∞
= 66 MPa and supershear rupture
velocity in the source region V
force
= 4330 m/s. The source excites
a strong wrinkle-like pulse in the positive direction with a velocity
€
V
r
+
≈ c
GR
and a weak pulse in the negative direction with a veloc-
ity
€
V
r
−
≈ α
−
.
30
2.14 Spatio-temporal evolution of slip velocity for a case similar to that
of Fig. 2.13 but with higher initial shear stress τ
∞
= 70 MPa.
31
3.1 Spatio-temporal evolution of slip velocity for a case similar to that
of Fig. 2.13 but with higher initial shear stress τ
∞
= 70 MPa.
40
ix
3.2 Effects of parameters of the employed rate- and state-dependent
friction. (a) The parameter V
0
controls the amplitude of the direct
response of the apparent coefficient of friction µ
app
to abrupt
changes of slip rate. Increasing V
0
reduces the direct effect. (b)
The parameter V
1
serves as a threshold slip rate for which µ
app
ap-
proaches the nominal dynamic level µ
d
. Increasing V
1
increases the
range of slip rates for which the frictional response is sensitive. (c)
The steady-state coefficient of friction as a function of slip rate for
different values of V
1
. (d) Evolution of µ
app
in response to abrupt
changes of slip rate under constant compressive normal stress for
two combinations of friction parameters used in the simulations.
45
3.3 Time histories of slip and shear stress (left columns) and func-
tional dependence of the shear stress on slip (right columns) at the
point x
1
= 1.0 m on the interface for (a) Case I with supershear
crack-like rupture, (b) Case II with subshear crack-like rupture, (c)
Case III with subshear single pulse, and (d) Case IV with super-
shear train of pulses.
52
3.4 Rupture front speeds V
R
in the positive x
1
direction versus time for
Cases I-IV. The speeds of the P wave (α), S wave (β) and √2b are
indicated by the dashed lines.
53
3.5 Profiles of frictional slip rate (red) and shear traction (blue) along
the interface at time 0.7 ms for (a) Case I, (b) Case II, (c) Case III,
and (d) Case IV. The insets show close-up views of the wiggles in
the profiles of slip rate.
55
3.6 Time histories of various energy components for (a) Case I, (b)
Case II, (c) Case III and (d) Case IV. The insets illustrate the ap-
proximate partition of the total change of elastic energy between
the kinetic energy W
kin
and heat energy W
heat
.
56
3.7 The mismatch in percentage between the energies
€
W
vol
el
+W
int
el
+W
heat
+W
kin
and
€
W
app
for Cases I-IV.
57
3.8 (a) The spatial distribution of shear stress σ
12
at time 0.5 ms and
(b) snapshots of the slip profiles along the interface at various
times for Case I with supershear crack-like rupture.
59
3.9 (a) The spatial distribution of shear stress σ
12
at time 0.7 ms and
(b) snapshots of the slip profiles along the interface at various
times for Case II with subshear crack-like rupture.
60
3.10 (a) The spatial distribution of shear stress σ
12
at time 1.5 ms and
(b) snapshots of the slip profiles along the interface at various
times for Case III with subshear pulse rupture.
61
x
3.11 (a) The spatial distribution of shear stress σ
12
at time 0.7 ms and
(b) snapshots of the slip profiles along the interface at various
times for Case IV with supershear train of pulses.
62
3.12 Simulations results with a reduced nominal dynamic coefficient of
friction µ
d
= 0.3 and the other parameters same as in Case II. (a)
Profiles of frictional slip rate (red) and shear traction (blue) along
interface at time 0.5 ms. (b) Rupture front speeds V
R
versus time.
(c) Snapshots of the slip profile along the interface at various
times. (d) Time histories of various energy components with inset
showing the approximate partition of the total change of elastic
energy between the kinetic energy W
kin
and heat energy W
heat
.
63
4.1 Geometry of the configuration analyzed and time function of the
impact loading. The dashed line denote interface S
int
(x
2
= 0)
where frictional sliding occur.
71
4.2 Close-up view of the mesh near the transition x
1
= 60 mm between
the region of fine uniform mesh and the region of graduated mesh
along S
int
.
73
4.3 Along-interface profiles of
€
Δ˙ u
slip
and T
s
/T
n
in Case I at 20, 25 and
30 µs for (a) the elastic calculation and (b) the plastic calculation.
79
4.4 Propagation speed of the frictional sliding front, V
slip
, versus time
in Case I for the elastic and plastic calculations.
81
4.5 Time histories of various energy components in Case I for (a) the
elastic calculation and (b) the plastic calculation.
82
4.6 Contours of Mises effective stress in Case I for (a) the elastic cal-
culation at t = 34.6 µs and (b) the plastic calculation at t = 34.6 µs.
The frictional sliding front is at x
1
= 58.0 mm in (a) and at x
1
=
55.5 mm in (b).
83
4.7 Along-interface profiles of
€
Δ˙ u
slip
and T
s
/T
n
in Case II at 30, 35 and
40 µs for (a) the elastic calculation and (b) the plastic calculation.
84
4.8 Propagation speed of the frictional sliding front, V
slip
, versus time
in Case II for the elastic and plastic calculations.
87
4.9 Time histories of various energy components in Case II for (a) the
elastic calculation and (b) the plastic calculation.
87
4.10 Contours of Mises effective stress in Case II for (a) the elastic cal-
culation at t = 37.6 µs and (b) the plastic calculation at t = 37.8 µs.
The frictional sliding front is at x
1
= 42.2 mm in (a) and at x
1
=
42.4 mm in (b).
88
xi
4.11 Along-interface profiles of
€
Δ˙ u
slip
and T
s
/T
n
in Case III (a) at 10, 20
and 29.4 µs for the elastic calculation and (b) at 10, 20, 30 µs for
the plastic calculation.
89
4.12 Propagation speed of the frictional sliding front, V
slip
, versus time
in Case III for the elastic and plastic calculations.
91
4.13 Time histories of various energy components in Case III for (a) the
elastic calculation and (b) the plastic calculation.
91
4.14 Contours of Mises effective stress in Case III for (a) elastic calcu-
lation at t = 29.4 µs and (b) the plastic calculation at t = 30 µs. The
frictional sliding front is at x
1
= 59.1 mm in (a) and at x
1
= 58.2
mm in (b).
92
4.15 Along-interface profiles of
€
Δ˙ u
slip
with the frictional slip front posi-
tion at 50 mm in (a) the elastic calculations and (b) the plastic cal-
culations. V
imp
= 20 m/s and C
s
= 100, 200 and 400 GPa/m. The
time for each profile is indicated.
93
4.16 Propagation speed of the frictional sliding front, V
slip
, versus time
for the elastic and plastic calculations with varying C
s
.
94
4.17 Along-interface profiles of
€
Δ˙ u
slip
and T
s
/T
n
at 10, 20 and 30 µs for
(a) the elastic calculation and (b) the plastic calculation with V
imp
=
5 m/s and all other parameters the same as for Case III.
96
4.18 Along-interface profiles of
€
Δ˙ u
slip
and T
s
/T
n
at 10, 20 and 28.5 µs
for (a) the elastic calculation and (b) the plastic calculation with
V
imp
= 30 m/s and all other parameters the same as for Case III.
96
4.19 Propagation speed of the frictional sliding front, V
slip
, versus time
for the elastic and plastic calculations with varying V
imp
.
97
4.20 The effect of impact velocity V
imp
on the time variations of (a) the
kinetic energy W
kin
and (b) the frictional dissipation along the in-
terface φ
slip
.
97
4.21 Contours of effective plastic strain for (a) the plastic calculation in
Case II with σ
0
= 27 MPa at t = 37.8 µs and (b) the plastic calcula-
tion in Case III with σ
0
= 27 MPa at t = 30 µs.
99
xii
4.22 Phase diagrams showing the effects of impact velocity V
imp
and
interface stiffness parameter C
s
on the mode of frictional sliding
with (a) elastic bulk material behavior and (b) plastic bulk material
behavior. Shaded symbols denote the combinations of parameters
that produce different frictional sliding modes in elastic and plastic
calculations. The initial compressive normal stress is Σ
0
= 10 MPa
and the reference flow strength for the plastic calculations is σ
0
=
10 MPa.
101
5.1 Cartesian and spherical-polar coordinate systems. The fault lies on
the x
1
-x
2
plane with the x
3
axis as the fault normal. For a shear dis-
location source the fault slips along the x
1
-axis direction and for a
tensile dislocation the fault displacement motion is along the x
3
-
axis direction.
108
5.2 Far-field radiation patterns of P-wave (top plots) and S-wave (bot-
tom plots) from (a) a shear dislocation source and (b) a tensile dis-
location source in a Poissonian solid (ν = 0.25). Fault-normal di-
rection and slip directions are labeled and ‘N’ denotes the neutral
axis.
109
5.3 (a) A model consisting of two half-spaces with the x-y plane as the
fault plane and the z-axis in the fault-normal direction. The dislo-
cation source (red star on the x-axis) can have both fault-parallel
(marked as 1) and fault-normal (marked as 2) components. x
s
, z
r
and R are the source depth with respect to the receiver plane (fixed
at 1 km), fault-normal distance and the hypocentral distance. (b)
Various receiver configurations on plane x = 0 employed in the
calculations of synthetic seismograms: Configuration ‘RCA’ with
receivers along a straight line y = 0 normal to the fault and through
the epicenter; Configuration ‘RCB’ with receivers along two
straight lines parallel (and very close) to the fault trace with z =
±0.1 m; Configuration ‘RCC’ with receivers along a circle cen-
tered at the epicenter with radius r = 0.2 km; and Configuration
‘RCD’ with receivers along a straight line y = 1 km. Receivers in
each configuration are in pairs symmetric about the fault plane.
110
5.4 Synthetic seismograms from a shear dislocation recorded with
RCD. The red and blue traces are obtained from numerical imple-
mentation of the analytical solution given by Ben-Zion [1990,
1999] whereas the black thin traces are obtained from numerical
implementation of the analytical solution given in Aki and Rich-
ards [2002].
111
xiii
5.5 The radiation curves (amplitude profiles of P and S phases) ob-
tained from theoretical formula (solid lines) and measurements
from synthetic seismograms (diamond and circle symbols) from a
shear dislocation in a homogeneous solid with (a) array RCA
where θ is defined in Fig. 5.2a, (b) array RCB and (c) array RCC
where α is the angle measured from +y axis clock-wise on the y-z
plane.
112
5.6 The radiation curves (amplitude profiles of P and S phases) that
are counterparts of plots in Fig. 5.5 from a tensile dislocation.
113
5.7 Three-component synthetic seismograms from a shear dislocation
in a homogeneous solid recorded with array (a) RCA (b) RCB and
(c) RCC. The vertical components in (a) are identically zero. Both
fault-parallel and vertical components in (b) are several orders of
magnitudes smaller than the fault-normal components. The verti-
cal components in (c) are amplified by 10 times. The largest wave-
form amplitudes are labeled on the right side in unit of m.
115
5.8 Three-component synthetic seismograms from a tensile dislocation
source in a homogeneous solid recorded with array (a) RCA (b)
RCB and (c) RCC. The vertical components in (a) and fault-
normal components in (b) are identically zero. The largest wave-
form amplitudes are labeled on the right side in unit of m.
116
5.9 The radiation curves (amplitude profiles of P and S phases) meas-
ured from synthetic seismograms for a shear dislocation between
dissimilar solids (10% velocity contrast) with array (a) RCA, (b)
RCB and (c) RCC. Solid and dashed curves denote measurements
on the fast-medium side and slow-medium side, respectively.
117
5.10 The radiation curves (amplitude profiles of P and S phases) meas-
ured from synthetic seismograms for a tensile dislocation between
dissimilar solids (10% velocity contrast) with array (a) RCA, (b)
RCB and (c) RCC. Solid and dashed curves denote measurements
on the fast-medium side and slow-medium side, respectively.
117
5.11 The spatial distribution of P- and S-wave amplitudes on a receiver
plane x = 0 radiated from a shear dislocation between (a) similar
and (b) dissimilar solids. The amplitudes are normalized with the
maximum value measured. The fault-normal distance (along z-
axis) and the along-strike distance (along y-axis) are both normal-
ized by source depth x
s
. The solid line at z = 0 denotes the fault
trace and the white dot denotes the epicenter. The head waves are
recorded as first arrivals in the region between the dashed line and
the solid line (fault).
119
xiv
5.12 The spatial distribution of amplitudes of P- and S-waves on a re-
ceiver plane x = 0 radiated from a tensile dislocation between (a)
similar and (b) dissimilar solids. The notations and symbols are
the same in Fig. 5.11.
120
5.13 Contours of the degree of expected interference between the head
and direct S waves based on analytical calculation of the differen-
tial travel-time between the different phases, source rise time, and
amplitude reduction due to geometrical spreading. The notations
and symbols are as in Fig. 5.11. See text for additional explana-
tion.
121
5.14 Three-component synthetic seismograms from a shear dislocation
between two dissimilar solids (10% velocity contrast) recorded
with array (a) RCA (b) RCB and (c) RCC. P and S head waves are
present on seismograms recorded within critical fault-normal dis-
tances. The vertical components in (c) are amplified by 10 times.
The largest waveform amplitudes are labeled on the right side in
unit of m.
124
5.15 Three-component synthetic seismograms from a tensile dislocation
between two dissimilar solids (10% velocity contrast) recorded
with array (a) RCA (b) RCB and (c) RCC. P and S head waves are
present on seismograms recorded within critical fault-normal dis-
tances. The vertical components in (c) are amplified by 10 times.
The largest waveform amplitudes are labeled on the right side in
unit of m.
125
xv
Abstract
This thesis examines dynamic ruptures along frictional interfaces and seismic radia-
tion in models of earthquake faults separating similar and dissimilar solids with the goal
of advancing the understanding of earthquake physics.
The dynamics of Mode-II rupture along an interface governed by slip-weakening fric-
tion between dissimilar solids are investigated. The results show that the wrinkle-like
rupture along such interfaces evolves to unilateral propagation in the slip direction of the
compliant side for a broad range of conditions, and the closer the initial shear stress is to
the static friction the smaller degree of material contrast is needed for this evolution to
occur. Transition of the wrinkle-like pulse to crack-like rupture occurs when the reduc-
tion of friction coefficient is sufficiently large.
Energy partition associated with various rupture modes along an interface governed
by rate- and state-dependent friction between identical solids is investigated. The rupture
mode changes with varying velocity dependence of friction, strength excess parameter
and length of the nucleation zone. High initial shear stress and weak velocity dependence
of friction favor crack-like ruptures, while the opposite conditions favor the pulse-like
mode. The rupture mode can switch from a subshear single pulse to a supershear train of
pulses when the width of the nucleation zone is increased. The elastic strain energy re-
leased over the same propagation distance by the different rupture modes has the order:
xvi
supershear crack, subshear crack, supershear train-of-pulses and subshear single-pulse.
General considerations and observations suggest that the subshear pulse and supershear
crack are, respectively, the most and least common modes of earthquake ruptures.
The effect of plasticity and interface elasticity on dynamic frictional sliding along an
interface induced by edge impact loading between two identical elastic-viscoplastic solids
is analyzed. The material on each side is isotropically strain-hardening. The interface is
characterized as having an elastic response together with an inelastic response character-
ized by rate- and state-dependent friction. The results show that bulk material plasticity
tends to smooth out oscillations. Larger impact velocity induces more extensive plastic
dissipation and larger effect on slip mode and energy partition. Also larger values of the
interface shear stiffness tend to favor crack-like mode of sliding.
The last part of we examine the characteristics of seismic radiation from localized
fault-opening and shear motions and the effect of having dissimilar solids across the fault
on seismic radiation by employing calculations of synthetic seismograms generated at
various receiver locations by shear and tensile dislocation sources. The existence of a ve-
locity contrast across the fault produces complexities that mask somewhat the polarity
and amplitude signals of body waves characteristic of the tensile dislocation. The record-
ing and analysis of the discussed signals for regular earthquakes that are dominated by
shear motion will require high-resolution receivers located very close to the fault.
1
Chapter 1
Introduction
Earthquakes are complex natural phenomena the physics of which is still not fully
understood despite significant progress that has been made in the past few decades. Fric-
tional sliding on fault plane and seismic radiation are two basic problems in earthquake
studies. This work uses numerical methods to examine the dynamics of frictional sliding
and the seismic radiation in models of earthquake faults separating similar and dissimilar
solids, with the goal of improving our understanding of earthquake physics.
The propagation of earthquake rupture is a large-scale frictional sliding process, the
dynamics of which can be highly nonlinear and complex with many small-scale features
playing significant roles. Traditional treatment of earthquake ruptures as cracks on planar
faults governed by near-constant friction between identical elastic solids under homoge-
neous stress loading is known to be oversimplified in many earthquake source studies
nowadays. Better understanding of earthquake source dynamics requires fault rupture
models with more realistic physical ingredients that have significant influences on earth-
quake rupture dynamics.
2
As shown by previous numerical and experimental studies, dynamic rupture along a
frictional interface can propagate in different modes, e.g., crack-like, pulse-like and
pulse-train, and at different speeds, e.g., subshear and supershear. Physical ingredients in
natural fault zones that can have significant effects on the properties of earthquake rup-
tures include: geometrical and constitutive properties of the fault, presence of fault bra-
ches, stress loading condition, elastic properties and plastic yielding or damage properties
of the rocks. A highly advanced fault model should further include the interplay among
all the above ingredients to account for the evolution of the fault, which at present is hard
to achieve. In this work, we incorporate some of those ingredients in our numerical mod-
els to investigate their effects on the dynamic rupture process.
In Chapter 2, we study the properties of the wrinkle-like pulses along an interface that
separates linear elastic solids with different elastic properties using finite difference
method. Theoretical studies show that mode-II ruptures along such frictional bimaterial
interface induce dynamic changes of normal stress (hereafter referred to as the ‘bimaterial
effect’). At subshear speed, rupture propagation in the direction of slip on the more com-
pliant side experiences dynamic reduction of normal stress leading to decrease of fault
strength, whereas rupture propagation in the other direction experiences dynamic increase
of normal stress leading to increase of fault strength. Previous numerical studies on wrin-
kle-like rupture assumed constant friction on the interface. However, the friction coeffi-
cient of rocks typically drops when slip occurs as observed in rock sliding experiments.
Here we employ a slip-weakening law to characterize the frictional property of the inter-
face. The reduction of fault strength due to the reduction of friction coefficient can add to
3
or lessen the bimaterial effect on the fault strength depending on the rupture propagation
direction along the bimaterial interface as described above. In our model the evolution
behavior of the wrinkle-like pulse to unilateral propagation with different degrees of ma-
terial contrast and closeness of the initial shear stress to static friction is examined. The
transition of rupture mode from wrinkle-like pulse to crack-like rupture with increasing
reduction of friction coefficient is observed. Rupture propagation speeds for pulse-like
and crack-like ruptures are also examined for different nucleation parameters.
During earthquake rupture, elastic strain energy stored in the crust by the tectonic
loading process is suddenly released and partitioned into other forms of energies includ-
ing seismic radiation energy and dissipative energies on the fault and in the bulk. Given
the lack of observational constraints on most of the energy components, a better under-
standing of the whole energy budget will benefit the development of an improved quanti-
tative theory for dynamic earthquake rupture and provide valuable information for seis-
mic hazard analysis related to strong motion. We study in Chapter 3 the energy partition
during dynamic in-plane frictional sliding along an interface governed by a rate- and
state-dependent friction between identical elastic solid using finite element method. In
our model ruptures propagate in different modes (crack-like, pulse-like or train-of-pulses)
and at different speeds (subshear or supershear) depending on the velocity dependence of
the friction, the strength excess parameter and the length of the nucleation zone. The en-
ergy partition pattern varies with different rupture properties that include supershear
crack-like rupture, subshear crack-like rupture, subshear single pulse and supershear train
of pulses. General considerations on the conditions that lead to the various rupture modes
4
together with observations suggest that the subshear pulse and supershear crack are, re-
spectively, the most and least common modes of earthquake ruptures.
In the dynamic rupture process, the elastically predicted stress concentration near the
rupture tip at the interface can be so high that the materials off the interface can yield ine-
lastically. By using the same finite element method as in Chapter 3 and a different model
configuration, we explore the effect of plasticity and interface elasticity on dynamic fric-
tional sliding between two identical elastic-viscoplastic solids in Chapter 4. The two sol-
ids are held together by a uniform compressive stress and dynamic frictional sliding is
initiated by shear loading by edge impact near the interface. As in Chapter 3, the interface
is characterized by an elastic response and an inelastic response characterized by a rate-
and state-dependent friction. The bulk materials are isotropically strain hardening. Effect
of plasticity is examined by comparing results from simulations with elastic material
properties and plastic material properties. Different sliding modes are obtained with dif-
ferent bulk material properties, interface elasticity and edge impact velocity. The associ-
ated energy partition patterns are also analyzed.
In Chapter 5, we examine the seismic radiation from localized fault-opening and
shear motion, along with the effect of having dissimilar solids across the fault on seismic
radiation. Earthquakes in geothermal, volcanic and other extensional regimes have tensile
components due to the volumetric changes in the source region. Also certain physical
processes during shear faulting can cause transient fault-opening motion, such as colli-
sion of geometric irregularities on rough surfaces, Schallamach-type waves and dynamic
change of normal stress due to the bimaterial effect as discussed in Chapter 2. In this last
5
chapter, synthetic seismograms from shear and tensile point-dislocation sources are ana-
lyzed and compared in search of seismic signatures characteristic of fault-opening mo-
tion. Synthetic seismograms from models of homogeneous whole space and dissimilar
half-spaces are also compared in order to investigate the effect of having different elastic
solids across the fault on seismic radiation.
6
Chapter 2
Dynamic Rupture Along a Bimaterial Interface
Governed by Slip-Weakening Friction
2.1 Introduction
Plate-boundaries and other major fault-zone structures have bimaterial interfaces that
separate rock with different mechanical properties. Such interfaces are created by large
motion across major faults and cumulative production of rock damage associated with the
faulting process [e.g., Ben-Zion and Sammis, 2003 and references therein]. For example,
seismic imaging studies based on body waves [e.g., Eberhart-Phillips and Michael, 1998;
Fuis et al., 2003; Thurber et al., 2004, 2006], surface waves [e.g., Shapiro et al., 2005]
and head waves [e.g., Ben-Zion and Malin, 1991; Ben-Zion et al., 1992; McGuire and
Ben-Zion, 2005; Zhao et al., 2009] as well as geodetic data modeling studies [e.g., Le
Pichon et al., 2005; Fialko et al., 2006; Wdowinski et al., 2007] indicate that the contrast
of seismic velocities at various locations across the San Andreas fault is about 5−30%.
Theoretical works indicate that mode II ruptures along a frictional bimaterial interface
between linear elastic solids produce dynamic changes of normal stress on the fault that
depend on the slip gradient, material properties and direction of rupture propagation [e.g.,
Weertman, 1980; Adams, 1995; Ben-Zion, 2001; Ranjith and Rice, 2001]. When sub-
7
sonic rupture propagates in the direction of slip on the compliant side of the fault (re-
ferred to below as the ‘positive’ direction), there is a dynamic reduction of normal stress
across the fault at the rupture tip, whereas for rupture propagation in the opposite direc-
tion (referred to as the ‘negative’ direction) there is a dynamic increase of normal stress
at the rupture tip. The magnitudes of these effects increase with the propagation velocity
and the degree of material contrast across the fault, up to the largest contrast for which
the generalized Rayleigh wave speed c
GR
exists. For typical contrasts of rock densities
and Poisson ratios, c
GR
exists up to S-wave velocity contrasts of about 30−40% and is
slightly below the slower S-wave speed [e.g., Weertman, 1963; Achenbach and Epstein,
1967; Ben-Zion, 2001; Ranjith and Rice, 2001]. In addition, the dynamic changes of
normal stress increase with propagation distance along the bimaterial interface due to a
dynamic instability [e.g., Adams, 1995, 1998; Ranjith and Rice, 2001] that produces a
continual transfer of energy to shorter wavelengths during rupture propagation [e.g., Ben-
Zion and Huang, 2002]. The observed contrasts of seismic velocities across major faults
fall in the range (e.g., 10−30%) for which the bimaterial effects are prominent.
Andrews and Ben-Zion [1997], Ben-Zion and Andrews [1998], Cochard and Rice
[2000] and Ben-Zion and Huang [2002] performed numerical simulations of ruptures
along a bimaterial interface governed by Coulomb and regularized Prakash-Clifton fric-
tion laws, with slip-independent friction coefficient, between two linear elastic solids.
The calculations have shown, in agreement with the above analytical works, that mode II
rupture along a bimaterial interface propagates for broad ranges of conditions in a self-
sustaining manner as a narrow ‘wrinkle-like’ pulse in the positive direction. In those stud-
8
ies the friction law is stable with respect to slip and frictional failure occurs only due to
dynamic changes of normal stress. This allows a clear focus on effects the assumed fault
zone structure has on dynamic rupture properties, without mixing effects of different fail-
ure mechanisms. In such cases, the rupture evolves quickly to a unilateral pulse that
propagates in the positive direction with a velocity close to c
GR
and a large dynamic re-
duction of normal stress at the propagating rupture tip. During propagation in the positive
direction the slip velocity increases and the width of the rupture pulse decreases, while
the amount of slip remains about constant along the fault.
Adams [2001], Ranjith and Rice [2001] and Cochard and Rice [2000] found that slip
pulses with a velocity near the slower P-wave speed α− can also propagate along a bima-
terial interface with a constant friction coefficient in the negative direction. However,
these pulses are more difficult to excite and considerably weaker than the primary wrin-
kle-like pulses in the positive direction. Brietzke and Ben-Zion [2005] simulated ruptures
in structures with multiple possible failure planes, two of which are bimaterial interfaces
on the opposite sides of a low-velocity fault zone layer, and found that ruptures tend to
migrate spontaneously to the bimaterial interfaces. In some cases slip pulses propagate
simultaneously in the two different positive directions on the opposite sides of the low-
velocity layer, producing together an effective bilateral rupture.
Anooshehpoor and Brune [1999] performed sliding experiments with two different
foam rubber blocks and found that large ruptures propagate only in the positive direction,
with properties compatible with the theoretical predictions for wrinkle-like slip pulses.
They also observed some small events with different propagation directions and various
9
transient properties. Xia et al. [2005] performed sliding experiments along a bimaterial
interface for several loading configurations and obtained asymmetric bilateral ruptures. In
all cases, the rupture fronts in the positive direction had stable properties that did not de-
pend on the employed experimental conditions and propagated with a velocity close to
c
GR
. In contrast, properties of ruptures in the negative direction varied from case to case,
and had transient features (e.g., transition from sub-Rayleigh to near α−) within given
experiments. McGuire et al. [2002] and Henry and Das [2001] analyzed rupture proper-
ties of over 100 large global earthquakes and found that most are predominantly unilat-
eral. High-resolution locations of several thousand small events on the San Andreas fault
show a directional asymmetry [Rubin and Gillard, 2000] that is compatible with a pre-
ferred propagation direction associated with the local velocity structure.
The above observational results suggest that earthquake ruptures in large fault zones
can have a preferred propagation direction that is correlated with the velocity structure as
predicted for wrinkle-like ruptures. As summarized by Ben-Zion [2001], the dynamic
properties of wrinkle-like rupture pulses can have important implications for effective
constitutive laws, expected ground motion, rise-time of earthquake slip, frictional heat,
suppression of branching in large fault zone structures, and various other issues of earth-
quake and fault mechanics. It is thus important to gain an improved understating of the
conditions leading to such ruptures and obtain additional related laboratory and field ob-
servations.
In this chapter, we attempt to clarify properties of mode II rupture along a bimaterial
interface governed by slip-weakening friction. We perform a systematic parameter-space
10
study associated with different values of the difference between static and dynamic fric-
tion coefficients, different values of initial shear stress, and different degrees of material
contrast across the fault. In the calculations with slip-weakening friction, the amount of
slip associated with the wrinkle-like pulses grows with propagation distance, in contrast
to the previous related results with slip-independent coefficient of friction. For differ-
ences between the static and dynamic coefficients of friction larger than about 0.5, rup-
ture along the bimaterial interface propagates as a bilateral crack with a velocity (in both
directions) close to c
GR
. Rupture pulses in the negative direction with velocity close to α−
are generated only when we use an initial imposed rupture velocity with a similar speed.
In such cases, the simulated rupture pulse in the positive direction still propagates with a
velocity close to c
GR
.
2.2 Analysis
2.2.1 Model setup
We consider 2D in-plane dynamic rupture problems along a frictional interface at
between two linear isotropic elastic half spaces (Fig. 2.1a). The calculations are carried
out using a second-order finite-difference formulation with frictional sliding on the fault
governed by the split-node sliding logic [e.g., Andrews, 1999]. The density, P- and S-
wave speeds of the two solids are ρ
i
, α
i
and β
i
, respectively, with subscripts i = 1 denoting
the material below the interface (y < 0) and i = 2 the material above (y > 0). The degree of
material contrast γ across the interface is defined as 1+γ = ρ
1
/ρ
2
= α
1
/α
2
= β
1
/β
2
. The ini-
tial stresses applied at the remote boundaries are
€
σ
xx
0
=σ
yy
0
=−σ
∞
and
€
σ
xy
0
=τ
∞
. With
11
these values, the initial maximum compressive principal stress direction is at 45º to the
fault plane.
Figure 2.1: (a) A Model for 2D finite-difference simulations of mode-II dynamic rupture
along a bimaterial interface. The initial normal and shear stress components at the remote
boundaries are −σ
∞
and τ
∞
. Events are nucleated by in two limited x−t domains that ex-
pand symmetrically about the origin. (b) A linear slip-weakening friction law that gov-
erns the propagation of spontaneous rupture outside the source region.
The simulations begin with initial shear stress along the fault that is below the static
frictional level everywhere along the fault. Ruptures are nucleated by gradually increas-
ing pore fluid pressure in two limited space-time domains that are symmetric about the
origin. This leads to a symmetrically expanding slip in a limited nucleation zone centered
at (x, y) = (0, 0). The boundaries of the imposed source are two ellipses (thick arrows near
the origin in Fig. 2.1a) in the positive and negative directions
€
1−ξ
2
−η
2
=0, where
€
ξ = x −vt
( )
a,
€
η = x +vt
( )
b−η
0
and
€
η
0
= a
2
+b
2
b. The pore pressure within these
two elliptical regions is given by
€
σ
p
= 1−ξ
2
−η
2
( )
σ
∞
, which rises smoothly from 0 at the
x
y
Nucleation zone
Material 2: !
2
, "
2
, #
2
Material 1: !
1
, "
1
, #
1
$%
&
&
$%
8
8
8
8
slip
friction coefficient
µ
d
s
D
c
µ
(a) (b)
12
ellipse boundaries to 1 at the ellipse centers. In all simulations we use a = 10 m, b = 60 m
and in most cases the rupture velocity in the imposed source region is v = V
force
= 2475
m/s (about 0.9 the Rayleigh wave speed of the stiffer material). In Section 2.2.6 a higher
value is used to facilitate the excitation of supershear ruptures.
Rupture propagation outside the imposed nucleation zone is governed by a slip-
weakening friction (Fig. 2.1b) with which the friction coefficient drops linearly from a
static value µ
s
to a dynamic level µ
d
over a characteristic slip distance D
c
. The size of the
spatial region behind the rupture tip where the amount of slip reaches D
c
is referred to as
the cohesive zone size X
c
. In dynamic rupture problems, X
c
shrinks due to a Lorentz-type
contraction as the rupture speed accelerates to a limiting value, which for mode-II rupture
is the generalized Rayleigh wave speed c
GR
[Weertman, 1980]. The values of material
contrasts employed in this study are within the range for which c
GR
exists. As mentioned
in the introduction, this range characterizes the contrast across natural faults. With a de-
gree of material contrast in this range and slip-weakening friction, the Adams instability
[Adams, 1995] exists for all values of the friction coefficient [Ranjith and Rice, 2001],
and the bimaterial effect increases with increasing degree of material contrast across the
fault [Ben-Zion and Andrews, 1998]. The existence of this Adams instability in our cal-
culations produces a transfer of energy to short wavelengths during rupture propagation,
which reduces dynamically the width of the rupture pulse [e.g., Ben-Zion and Huang,
2002].
Proper numerical calculations require the usage of a grid size dx that is small enough
to resolve the smallest existing physical length scale. For the problem at hand with evolv-
13
ing length scales, this is highly challenging and can be done at present only over a limited
propagation distance. Based on preliminary numerical calculations, a grid size dx = 0.5 m
was found to provide stable numerical results over a propagation range of 600 m for
many sets of conditions. In some cases, the calculations develop strong numerical oscilla-
tions over that propagation distance and can be interpreted only tentatively. The numeri-
cal time step is
€
dt =0.48 3dx α
+
with α
+
being the faster P-wave speed. In all calcula-
tions, α
+
= 5196.0 m/s and the numerical time step is 80 µs. Table 2.1 summarizes the
parameters that are common to all the cases examined in the paper.
As discussed by Brietzke and Ben-Zion [2005], the inclusion of artificial viscosity
can reduce the numerical oscillations [Harris and Day, 1997], but it may also suppress
genuine small-scale features and change additional aspects of the problem by absorbing
energy. To avoid altering the physics in this problem, we do not include artificial viscos-
ity in our simulations.
Parameter Value
Poisson’s ratios ν
1
, ν
2
0.25
Density of the stiffer material ρ
1
3333.3 kg/m
3
Faster S-wave speed 3000 m/s
Faster P-wave speed 5195 m/s
Static coefficient of friction µ
s
0.75
Characteristic slip distance D
c
0.1 m
Initial compressive normal stress σ
∞
100 MPa
Source parameter a 10 m
Source parameter b 60 m
Grid size dx 0.5 m
Time step dt 80 µs
Table 2.1: Values of model parameters common to all calculations.
14
2.2.2 A homogeneous case
Figure 2.2: (a) Profiles of slip (top panel) and slip velocity (bottom panel) along the fault
at different times for a homogeneous case. (b) The corresponding spatio-temporal evolu-
tion of slip velocity. The results are symmetric in the positive and negative directions.
To provide a reference for comparison with results of rupture along a bimaterial inter-
face, we first study properties of rupture in a homogeneous solid (γ = 0). The initial shear
stress is τ
∞
= 66 MPa and the dynamic coefficient of friction is µ
d
= 0.65. Fig. 2.2(a) gives
the calculated slip (top panel) and slip velocity (bottom panel) at different times and Fig.
0.02
0.04
Slip (m)
Homogeneous case
50dt
100dt
200dt
400dt
800dt
1200dt
1600dt
!400 !300 !200 !100 0 100 200 300 400
0
5
10
Position along fault (m)
Slip velocity (m/s)
0
V r = 2750 m/s V r = 2750 m/s
(slip velocity, m/s)
!400 !300 !200 !100 0 100 200 300 400
Position along fault (m)
0
0.02
0.04
0.06
0.08
0.1
0.12
0.14
Time (s)
8 9 7 6 5 4 3 2 1
! !
(a)
(b)
15
2.2(b) shows the spatio-temporal evolution of the slip velocity along the fault. In this ho-
mogeneous case, there is no coupling between the slip and changes of normal stress on
the fault and the rupture propagates, as expected, in a symmetric bilateral fashion. The
dark blue regions in Fig. 2.2(b) are associated with the imposed source nucleation and
weak rupture fronts generated by the radiated seismic waves. The main rupture fronts
propagate in both directions with a velocity close to the generalized Rayleigh wave speed
until they are arrested at ±3.8S from the origin, where S ~ 60 m is the source size. In-
creasing the initial shear stress τ
∞
and/or using larger differences of µ
s
– µ
d
produce simi-
lar bilateral ruptures that propagate larger distances along the fault.
2.2.3 Corresponding cases with material contrast
When the fault is a bimaterial interface, spatial variations of slip produce changes of
normal stress within the slip region [Weertman, 1980]. For subsonic propagation, this
creates dynamic dilatation at the rupture tip that propagates in the direction of slip on the
compliant side of the fault, which is the positive x direction in our coordinate system, and
dynamic compression at the rupture tip in the opposite direction. In this section we exam-
ine effects of material contrasts across the fault by using the same values of initial shear
stress τ
∞
= 66 MPa and dynamic coefficient of friction µ
d
= 0.65 as for the previous ho-
mogeneous case. In the following sections we consider variable values of τ
∞
and µ
d
.
Fig. 2.3(a) gives the slip and slip velocity along the fault for a 4% material contrast
across the interface. In this case the symmetrically initiated bilateral rupture evolves after
some propagation distance along the bimaterial interface to a unilateral rupture in the
positive direction. The corresponding spatio-temporal evolution of slip velocity is plotted
16
in Fig. 2.3(b). Rupture in the positive direction propagates at a speed close to c
GR
, while
rupture in the negative direction propagates at a slightly lower velocity. With the em-
ployed small γ = 0.04, µ
s
– µ
d
= 0.1 and initial stress configuration, the rupture dies in
both directions. However, the rupture in the positive direction has larger propagation dis-
tance, larger slip and larger slip velocity.
Figure 2.3: (a) Profiles of slip (top panel) and slip velocity (bottom panel) along the fault
at different times for a case with 4% material contrast. (b) The corresponding spatio-
temporal evolution of the slip velocity.
(a)
(b)
0.02
0.04
50dt
100dt
200dt
400dt
800dt
1200dt
1600dt
!400 !300
!200
!100 0 100 200 300 400
0
5
10
Position along fault (m)
Slip velocity (m/s)
0
Slip (m)
4% contrast
= 2700 m/s
(slip velocity, m/s)
!400 !300 !200 !100 0 100 200 300 400
Position along fault (m)
0
0.02
0.04
0.06
0.08
0.1
0.12
0.14
Time (s)
8 9 7 6 5 4 3 2 1 10
V r
!
V r = 2690 m/s
!
17
When the degree of material contrast is increased to 5%, the strength of the bimaterial
effects increases, and the rupture asymmetry becomes more pronounced (Fig. 2.4). In this
case the rupture dies more quickly in the negative direction, while it continues to grow in
the positive direction. With the employed frictional and stress parameters, the bimaterial
effects with γ = 0.05 are large enough to produce a self-sustaining unilateral rupture pulse
in the positive direction. The rupture velocities in the two opposite directions for this and
several other representative cases are listed in Table 2.2. Fig. 2.5 shows the evolution of
the slip profiles along the fault for cases with 10% contrast (top panel) and 20% contrast
(bottom panel). The asymmetry in the rupture properties for such cases is highly pro-
nounced. Note here the amount of slip increases with propagation distance in the positive
direction, in contrast to the results with a constant coefficient of friction [Andrews and
Ben-Zion, 1997; Ben-Zion and Andrews, 1998].
Figure 2.4: Profiles of slip and slip velocity along the fault at different times for a case
with 5% material contrast across the interface. The rupture is arrested quickly in the
negative direction and becomes self-sustained in the positive direction.
400 300 200 100 0 100 200 300 400
0.02
0.04
50dt
100dt
200dt
400dt
800dt
1200dt
1600dt
!400 !300 !200 !100 0 100 200 300 400
0
10
20
Position along fault (m)
Slip velocity (m/s)
0
Slip (m)
5% contrast
18
Model parameters Rupture velocity
γ
τ
∞
(MPa)
µ
d
V
force
(m/s)
c
GR
(m/s)
Figure #
€
V
r
+
(m/s)
€
V
r
−
(m/s)
0 66 0.65 2475 2758 2 2750 2750
0.04 66 0.65 2475 2703 3 2700 2690
0.05 66 0.65 2475 2689 2700 2625
0.1 66 0.65 2475 2620 2625 2560
0.2 66 0.65 2475 2475 2500 2490
0.2 66 0.05 2475 2475 11 5100
*
4300
*
0.2 66 0.75 3500 2475 2470 4230
0.2 66 0.60 4330 2475 2500 4240
0.2 66 0.65 4330 2475 13 2460 4250
0.2 66 0.75 4330 2475 2450 4250
0.2 68 0.65 4330 2475 2480 4300
0.2 68 0.55 5000 2475 2470 4310
0.2 70 0.65 4330 2475 14 2475 4330
Table 2.2: Rupture velocities for cases with different parameters. The symbol ‘*’ denotes
limit values of the rupture speed. Figure numbers denote plots showing rupture velocities
for that case.
Figure 2.5: Profiles of slip along the fault at different times for cases with 10% (top) and
20% (bottom) material contrast across the interface. The ruptures propagate in the posi-
tive direction in a self-sustaining manner with increasing slip.
400 300 200 100 0 100 200 300 400
0.05
0.1
50dt
100dt
200dt
400dt
800dt
1200dt
1600dt
!400 !300 !200 !100 0 100 200 300 400
0
0.1
0.2
Position along fault (m)
Slip (m)
20% contrast
0
Slip (m)
10% contrast
19
Figure 2.6: Maximum propagation distances of ruptures in the opposite directions vs. the
degree of material contrast. Cases with self-sustained rupture in the positive direction are
marked with “Divergence”. The inset gives the ratios of the maximum propagation dis-
tance of rupture in the positive direction divided by the corresponding distance in the
negative direction.
One way of quantifying the degree of asymmetry of rupture properties is to show the
maximum propagation distance in the two opposite directions for different sets of pa-
rameters. This is done in Fig. 2.6 with units of the imposed source size (S ~ 60 m) for τ
∞
= 66 MPa, µ
d
= 0.65, and 0 ≤ γ ≤ 0.2. Cases in which the rupture in the positive direction
becomes self-sustained are marked with the label “Divergence”. As γ increases, the rup-
ture in the negative direction is arrested earlier, while it propagates farther in the positive
direction and becomes self-sustained when the degree of contrast is between 4% and 5%.
The ratios of the maximum propagation distance in the positive direction divided by the
corresponding distance in the negative direction are given in the inset. The curve in the
0 5 10 15 20
0
1
2
3
4
5
6
7
8
Material contrast (%)
Final propagation distance / source size
negative dirction
positive direction
Divergence
0 1 2 3 4 5
1
1.5
2
2.5
3
3.5
4
Material contrast (%)
Ratio D+/D!
Divergence
8
8
µ
s
! µ
d
= 0.1, µ
s
" ! # = 9 MPa
20
inset shows clearly that the degree of rupture asymmetry grows rapidly with increasing
contrast across the fault.
A more physical way of quantifying the degree of rupture asymmetry is to plot for
various cases the scalar potency values [e.g., Ben-Zion, 2003], given by the integral of
slip over propagation distance, in both directions. This is done in the following sections,
where we show results for cases with different values of dynamic coefficient of friction
and initial shear stress.
2.2.4 Variable difference between static friction and initial shear stress
As the initial shear stress τ
∞
approaches the value of the static friction µ
s
σ
∞
along the
fault, the rupture can propagate, once initiated, for larger distances along the fault and
with higher slip velocity. In the previous two sections, the difference µ
s
σ
∞
− τ
∞
was 9
MPa. Here we keep the dynamic coefficient of friction µ
d
= 0.65 as in the previous sec-
tions, while varying the difference µ
s
σ
∞
− τ
∞
in the range 5−15 MPa.
Fig. 2.7 shows the scalar potency values generated by ruptures in the positive (red
lines) and negative (blue lines) directions for cases with different degrees of material con-
trast and initial shear stress values of τ
∞
= 62, 65 and 68 MPa. The corresponding ratios
of potencies in the positive and negative directions are given in the inset. The short verti-
cal dashed lines indicate the approximate values of material contrast beyond which the
final rupture propagation distance and potency value become unbounded in the positive
direction. The overall trends of potency values associated with the different cases are
similar to the evolution of maximum propagation distances shown in Fig. 2.6. With
higher initial shear stress, the rupture propagates for a larger distance in the negative di-
21
rection before being arrested, and it becomes self-sustained (associated with unbounded
propagation distance and potency) for a smaller degree of material contrast. For the cases
examined in Fig. 2.7, the rupture velocity is always close to c
GR
in the positive direction
and somewhat slower in the negative direction (Table 2.2). In addition, simulations with
larger initial shear stress produce more high-frequency numerical oscillations. This im-
plies (as might be expected) that the Adams instability [Adams, 1995] becomes more
vigorous as τ
∞
approaches the static friction level. Due to this problem, we cannot obtain
reliable results with values of τ
∞
greater than 70 MPa, i.e., µ
s
σ
∞
− τ
∞
< 5 MPa.
Figure 2.7: Scalar potencies of ruptures in the opposite directions versus degrees of mate-
rial contrast for several different values of initial shear stress. The inset shows the ratios
of rupture potency in the positive direction divided by the corresponding value in the
negative directions.
0 5 10 15 20
0.1
0.15
0.2
0.25
0.3
0.35
0.4
0.45
0.5
Material contrast (%)
Rupture potency (kmx cm)
= 68 MPa
= 65 MPa
= 62 MPa
!
0 1 2 3 4 5 6 7 8
1
1.1
1.2
1.3
1.4
1.5
1.6
1.7
1.8
1.9
2
Material contrast (%)
Ratio P+/P"
! = 68 MPa
! = 65 MPa
! = 62 MPa
µ
s
" µ
d
= 0.1
8
!
8
!
8
8 8 8
22
Figure 2.8: A phase diagram illustrating the effects of µ
s
σ
∞
− τ
∞
and degree of material
contrast on symmetry properties of ruptures along a bimaterial interface. The limit values
of 0 (deep blue) and 2 (dark red) represent, respectively, symmetric bilateral and fully
asymmetric unilateral ruptures. Detailed results associated with parameter values marked
by the white dots on the black horizontal line are shown in figures with numbers indi-
cated in the circles. The label ‘5t’ stands for the top panel of Fig. 2.5.
Shown in Fig. 2.8 is a phase diagram that summarizes symmetry properties of rup-
tures for cases with µ
s
– µ
d
= 0.1, µ
s
σ
∞
− τ
∞
in the range 5−15 MPa and different degrees
of material contrast. The color denotes the ratio of potency values in the positive direc-
tion divided by corresponding values in the negative direction. The color scale is chosen
to highlight details of the transition of ruptures from bilateral to predominately unilateral
propagation. Toward this end, potency ratios that are greater than or equal to 2 are treated
as 2−ε, where ε is a small positive number, while the value of 2 (dark red) is reserved for
cases in which the potency ratio diverges and the rupture becomes fully unilateral. A po-
23
tency ratio of one (deep blue) implies symmetric bilateral rupture over the range of calcu-
lated propagation distance (up to 600 m). A horizontal extension of the phase diagram to
the right for 0.1 < γ < 0.2 is associated entirely with values of 2 (unilateral ruptures) and
is not shown.
When the difference between the static strength and initial shear stress is 5 MPa, the
rupture evolves to a unilateral propagation with µ
s
– µ
d
= 0.1 for a material contrast that is
as small as about 2%. When the material contrast is about 10% or larger (up to the value
for which c
GR
exists), the rupture becomes unilateral for all examined cases of µ
s
σ
∞
− τ
∞
.
The transition from bilateral to unilateral rupture for different degrees of material contrast
occurs at a zone that scales approximately linearly (over the range of calculated propaga-
tion distance) with the difference µ
s
σ
∞
− τ
∞
.
2.2.5 Variable difference between the static and dynamic coefficients of friction
As shown in the previous sections, increasing degree of material contrast and increas-
ing initial shear stress promote the development of self-sustaining wrinkle-like ruptures
in the positive direction. Lager material contrast produces stronger dynamic reduction of
normal stress at the propagating tip in the positive direction, while higher initial shear
stress makes it easier for the rupture to continue to induce failure ahead. With a given
source properties and maximum propagation distance, the overall energetics of the rup-
ture is controlled by the difference between the static and dynamic coefficients of fric-
tion. Increasing this difference produces a larger static stress drop, resulting in larger slip
and potency values per rupture length. A large reduction of the friction coefficient, which
occurs in both propagation directions, may overcome the dynamic increase of normal
24
stress in the negative direction, and sustain the occurrence of bilateral ruptures over a
given propagation distance and other conditions.
Figure 2.9: A phase diagram illustrating the effects of µ
s
– µ
d
and degree of material con-
trast on symmetry properties of ruptures along a bimaterial interface. The other symbols
and notations are the same as in Fig. 2.8.
For the calculations in this section, the initial shear stress is fixed (as in Sections 2.2.2
and 2.2.3) at τ
∞
= 66 MPa, while the dynamic coefficient of friction varies from 0.75 to 0
(or equivalently the difference µ
s
– µ
d
varies from 0 to 0.75). Fig. 2.9 shows a phase dia-
gram of symmetry properties of ruptures for cases of µ
s
– µ
d
ranging from 0 to 0.35. The
color scale and other symbols are similar to Fig. 2.8. The results show a transition to uni-
lateral ruptures at about 5% contrast for all values of µ
s
– µ
d
in this examined range. Fig.
2.10 shows slip velocities on the fault at different times for 20% material contrast and
representative values of µ
s
– µ
d
in the range 0 to 0.75. For small values of µ
s
– µ
d
, the re-
25
sults are stable (over the calculated distance range) and show rapid rupture evolution to a
unilateral wrinkle-like pulse in the positive direction. As the difference µ
s
– µ
d
increases
toward 0.4 (Fig. 2.10c), high-frequency oscillations become very prominent and the qual-
ity of the results degrades. The high-frequency oscillations reflect the presence of nu-
merical noise, but may also be viewed as increasing sensitivity to a physical instability
[e.g., Adams, 1995; Ben-Zion, 2001] and hence approaching transition to a different dy-
namic regime. The oscillations may be interpreted as a train of pulses in an evolutionary
stage between dynamic regimes associated with pulse and crack modes of rupture [e.g.,
Coker et al., 2005]. This is supported by Figs 2.10(e)-(f), where larger values of µ
s
– µ
d
are seen to pro-duce smaller amounts of high-frequency oscillations, and the slip veloci-
ties exhibit a transition from unilateral wrinkle-like pulses to bilateral crack-like ruptures.
We note that along with the transition in the symmetry properties and mode of ruptures,
the maximum slip velocity decreases from Fig. 2.10(a) to Fig. 2.10(f).
The slip velocity in Fig. 2.10(f) with µ
s
– µ
d
= 0.75 is larger in the negative direction
than in the positive direction. This suggests that the rupture velocity in this case is super-
shear, in which case the bimaterial effects are reversed and there is reduction of normal
stress in the negative direction and increase in the positive direction [Weertman, 2002].
The spatio-temporal evolution of slip velocity for a bilateral crack-like rupture generated
with µ
s
– µ
d
= 0.7, τ
∞
= 66 MPa, and γ = 0.2 is shown in Fig. 2.11. In this case, the rupture
accelerates rapidly to supershear speeds in both directions and then gradually approaches
the slower P-wave speed in the negative direction α
−
and the faster P-wave speed in the
26
Figure 2.10: Profiles of slip velocity along the fault at different times for 20% material
contrast across the interface and different values of µ
s
– µ
d
. The high-frequency oscilla-
tions near µ
s
– µ
d
= 0.4 are associated with a transition from unilateral wrinkle-like rup-
ture pulses (a-b) to bilateral crack-like ruptures (e-f).
!400 !300 !200 !100 0 100 200 300 400
0
20
40
60
80
100
120
Position along fault (m)
Slip velocity (m/s)
50dt
100dt
200dt
300dt
400dt
800dt
1000dt
1200dt
!400 !300 !200 !100 0 100 200 300 400
0
10
20
30
40
50
60
70
80
90
100
Position along fault (m)
Slip velocity (m/s)
50dt
100dt
200dt
300dt
400dt
800dt
1000dt
1200dt
!400 !300 !200 !100 0 100 200 300 400
0
20
40
60
80
100
120
Slip velocity (m/s)
50dt
100dt
200dt
300dt
400dt
800dt
1000dt
1200dt
Position along fault (m)
(a)
(b)
(c)
µ
s
= 0.75, µ
d
= 0.70
µ
s
= 0.75, µ
d
= 0.50
µ
s
= 0.75, µ
d
= 0.35
27
Figure 2.10: (Continued.)
(d)
(e)
(f)
!400 !300 !200 !100 0 100 200 300 400
0
10
20
30
40
50
60
70
80
Slip velocity (m/s)
50dt
100dt
200dt
300dt
400dt
800dt
1000dt
1200dt
Position along fault (m)
!400 !300 !200 !100 0 100 200 300 400
0
5
10
15
20
25
30
35
40
45
Position along fault (m)
Slip velocity (m/s)
50dt
100dt
200dt
300dt
400dt
800dt
1000dt
1200dt
!400 !300 !200 !100 0 100 200 300 400
0
5
10
15
20
25
30
35
40
45
Position along fault (m)
Slip velocity (m/s)
50dt
100dt
200dt
300dt
400dt
800dt
1000dt
1200dt
µ
s
= 0.75, µ
d
= 0.25
µ
s
= 0.75, µ
d
= 0.10
µ
s
= 0.75, µ
d
= 0.
28
Figure 2.11: Spatio-temporal evolution of slip velocity for a case with µ
s
– µ
d
= 0.7. Su-
pershear crack fronts propagate in both directions with velocities close to the P-wave
speed. The crack front in the positive direction is followed by a slip pulse with a velocity
close to c
GR
.
positive direction α
+
. The initial crack rupture front in the positive direction is followed
by a wrinkle-like pulse with a velocity
€
V
r
+
≈ c
GR
.
Fig. 2.12 summarizes the symmetry properties of ruptures for different values of µ
s
–
µ
d
, different degrees of material contrast and µ
s
σ
∞
− τ
∞
= 9 MPa. The range µ
s
– µ
d
be-
tween 0.35 and 0.65 is left blank since the high-frequency oscillations in this range pre-
vent us from obtaining reliable potency values. The dashed line in the phase diagram in-
dicates a tentative transition that separates dynamic regimes associated with predomi-
nately bilateral and predominately unilateral raptures. Cases with µ
s
– µ
d
< 0.35 lead to
narrow pulses (self-sustained in the positive direction and decaying in the negative one)
= 4100 m/s
= 4300 m/s
= 3750 m/s
= 5100 m/s
V = 2470 m/s r
(slip velocity, m/s)
!400 !300 !200 !100 0 100 200 300 400
Position along fault (m)
0
0.02
0.04
0.06
0.08
0.1
0.12
0.14
Time (s)
5 10 15 20 25 30 35 40
Vr
!
Vr
!
Vr
!
Vr
!
!
29
that propagate with
€
V
r
+
≈ c
GR
, while cases with µ
s
– µ
d
> 0.6 that lead to bilateral crack-
like ruptures have supershear speeds similar to the results of Fig. 2.11.
Figure 2.12: A summary diagram extended from Fig. 2.9. The dashed line indicates a ten-
tative transition that separates dynamic regimes associated with predominantly bilateral
and predominantly unilateral ruptures.
2.2.6 Supershear ruptures with a larger V
force
As indicated in the previous sections, our simulations with a gradual imposed source
that expands with subshear rupture velocity, and are associated with µ
s
– µ
d
< 0.4 that
leads to unilateral wrinkle-like ruptures, do not produce supershear rupture in either
propagation direction (Table 2.2). A supershear rupture may be induced by using a more
vigorous source, which can be achieved with our nucleation procedure by increasing the
source size or increasing the rupture velocity of the imposed source. In this section we
30
present simulations that use an imposed source with same size as before (Table 2.1) but
with an expanding velocity of V
force
= α
−
= 4330 m/s.
Figure 2.13: Spatio-temporal evolution of slip velocity for a case with 20% contrast, ini-
tial shear stress τ
∞
= 66 MPa and supershear rupture velocity in the source region V
force
=
4330 m/s. The source excites a strong wrinkle-like pulse in the positive direction with a
velocity
€
V
r
+
≈ c
GR
and a weak pulse in the negative direction with a velocity
€
V
r
−
≈ α
−
.
Fig. 2.13 shows the spatio-temporal evolution of the slip velocity for a case with 20%
material contrast across the interface, τ
∞
= 66 MPa and µ
d
= 0.65. In this and other cases
with small µ
s
– µ
d
and high V
force
, there is a stronger rupture pulse that propagates in the
positive direction with
€
V
r
+
≈ c
GR
and a weaker pulse that propagates in the negative direc-
tion with
€
V
r
−
≈ α
−
. If we increase the initial shear stress to 70 MPa while keeping V
force
and the other parameters the same, we obtain a pulse in the negative direction that propa-
gates with
€
V
r
−
≈ α
−
over the entire calculation range (Fig. 2.14). However, the amplitudes
of slip and slip velocity in the negative direction are again considerably smaller than
= 2460 m/s
= 4250 m/s
(m/s)
(slip velocity, m/s)
!400 !300 !200 !100 0 100 200 300 400
Position along fault (m)
0
0.02
0.04
0.06
0.08
0.1
0.12
0.14
Time (s)
10 20 30 40 60 70 80 50
Vr
!
Vr
!
31
those associated with the primary wrinkle-like pulse in the positive direction. Simulations
with values of V
force
between 3500 m/s and 5000 m/s (≈ α
+
) produce similar results (Ta-
ble 2.2).
Figure 2.14: Spatio-temporal evolution of slip velocity for a case similar to that of Fig.
2.13 but with higher initial shear stress τ
∞
= 70 MPa.
2.3 Discussion
We performed a detailed parameter-space study aiming to clarify properties of dy-
namic rupture on a bimaterial interface governed by slip-weakening friction. The study
focuses on effects of the (1) amount of reduction of friction coefficient with slip, (2) de-
gree of material contrast across the interface, and (3) closeness of the initial shear stress
to the static friction. Using 2-D finite difference calculations, we construct phase dia-
grams that show the effect of the above three factors on the transition of bilateral to uni-
lateral rupture propagation and on the transition of pulse-like to crack-like rupture. All
(slip velocity, m/s)
!400 !300 !200 !100 0 100 200 300 400
Position along fault (m)
0
0.02
0.04
0.06
0.08
0.1
0.12
0.14
Time (s)
20 40 60 80 120 140 100
= 2475 m/s
= 4330 m/s Vr
!
Vr
!
32
simulations begin with an imposed bilateral rupture in a limited source region that ex-
pands symmetrically with a velocity V
force
. In most cases, V
force
= 2475 m/s = c
GR
for 20%
contrast of shear wave velocities across the fault. This value is also about 0.9 the
Rayleigh wave speed of the stiffer material and ~0.9c
GR
for a few percent contrast across
the fault (e.g., γ < 0.05). A propagation speed of c
GR
is a stable rupture velocity along a
bimaterial interface [e.g., Weertman, 1963, 1980] and is at the high end of rupture veloci-
ties that are inferred for most large earthquakes [e.g., Ben-Zion, 2003]. The employed
smoothly expanding source is constructed [Andrews and Ben-Zion, 1997] to gradually
approach conditions that may be typical to the ensuing spontaneous ruptures.
The bimaterial effects (dynamic changes of normal stress at the rupture tips) increase
with propagation distance along the bimaterial interface. It is thus important to perform
calculations over large propagation distances that would allow separation of early tran-
sient effects from later (ideally limit) evolutionary behavior. However, the dynamic re-
duction of physical length scales in this problem, associated with Lorentz-like contraction
and the Adams instability, produces with certain parameters strong numerical oscillations
(Fig. 2.10), and limits the propagation distance (or equivalently number of calculation
steps) and parameter values that can be simulated reliably with the employed code [An-
drews and Ben-Zion, 1997; Andrews, 1999; Ben-Zion and Huang, 2002].
Within current computational limitations, we find that the wrinkle-like pulse has a
preferred propagation direction that is the same as the slip direction on the compliant side
of the fault. This result is compatible with expectations based on analytical studies
[Weertman, 1980; Adams, 1995] and previous numerical calculations employing slip-
33
independent friction coefficient [Andrews and Ben-Zion, 1997; Ben-Zion and Andrews,
1998; Cochard and Rice, 2000; Ben-Zion and Huang, 2002; Ben-Zion and Shi, 2005].
However, in contrast to previous results associated with those calculations, the incorpora-
tion of slip-weakening friction leads to slip values that increase with propagation distance
along the fault.
With the employed nucleation procedure, the degree of material contrast needed to
produce a transition from bilateral to unilateral wrinkle-like rupture does not depend on
the value of µ
s
– µ
d
(Fig. 2.9) in the range 0−0.3. However, if µ
s
– µ
d
is larger than about
0.5, there is a transition from wrinkle-like slip pulses to bilateral crack-like ruptures (Fig.
2.10). The transition near µ
s
– µ
d
≈ 0.5 is associated with high-frequency oscillations that
might be interpreted as a train of pulses. As the difference µ
s
– µ
d
increases further, the
amount of high-frequency oscillations decreases and the ruptures propagate over the cal-
culated distance as stable bilateral cracks (Figs 2.10-2.12).
As mentioned in the introduction, bimaterial interfaces with degrees of contrast in the
range for which c
GR
exists are present in the structures of large faults, and within this
range the bimaterial effects increase with the degree of material contrast. For given initial
shear stress and friction coefficients, the simulated wrinkle-like ruptures become more
energetic, with higher slip and slip velocity, as the degree of material contrast increases
(up to the value for which c
GR
exists). As might be expected, the ruptures also become
more energetic as the initial shear stress τ
∞
is increased toward the static friction f
s
σ
∞
. The
simulations show further that in cases leading to wrinkle-like pulses, the transition of rup-
34
ture from bilateral to unilateral propagation occurs (Fig. 2.8) with smaller difference µ
s
σ
∞
− τ
∞
for lower degrees of material contrast γ.
In cases for which the bimaterial effects are significant and the wrinkle-like rupture
evolves to a unilateral propagation, the rupture velocity in the positive direction
€
V
r
+
is
always near c
GR
(Table 2.2). In such cases, the propagation velocity in the negative direc-
tion
€
V
r
−
(before the pulse dies out) is also very close to c
GR
, as long as V
force
in the im-
posed source region is close to c
GR
. Simulations with imposed rupture speed in the source
region close to the α
−
wave velocity can excite a rupture pulse in the negative direction
with
€
V
r
−
≈ α
−
. These results are in agreement with previous numerical calculations [Har-
ris and Day, 1997; Cochard and Rice, 2000] and analytical results [Adams, 2001; Ranjith
and Rice, 2001]. When the parameters are such that the bimaterial effects are not signifi-
cant (and the initial shear stress is close to the static friction), we obtain supershear rup-
tures with velocities between the fast and slow P-wave speeds. These results are com-
patible with previous theoretical and observational results for rupture in a homogeneous
solid [e.g., Andrews, 1976; Rosakis et al., 1999; Abraham and Gao, 2000; Ravi-Chandar
et al., 2000; Kubair et al., 2002; Coker et al., 2005].
Our parameter-space study provides a framework that can be used to understand the
range of results obtained in previous numerical simulations [Andrews and Ben-Zion,
1997; Ben-Zion and Andrews, 1998; Harris and Day, 1997; Cochard and Rice, 2000;
Ben-Zion and Huang, 2002] and related laboratory experiments of sliding along bimate-
rial interfaces [Anooshehpoor and Brune, 1999; Xia et al., 2005]. In the experiment of
Anooshehpoor and Brune [1999], the frictional and material properties were apparently in
35
the range for which the bimaterial effects are significant, and the gradual loading pro-
duced wrinkle-like ruptures that propagated in the positive direction with a velocity
€
V
r
+
≈
c
GR
. In the experiment of Xia et al. [2005], the bimaterial effects were apparently less
pronounced and the abrupt loading in the imposed source excited in some cases a pulse
that propagated in the negative direction with
€
V
r
−
≈ α
−
. Laboratory data that target broader
ranges of parameters will provide additional important constraints. For earthquake appli-
cations, it would be useful to use an initiation procedure that approximates the conditions
leading to typical large earthquake ruptures (i.e., gradual loading or initial imposed rup-
ture with subshear velocity). In addition, it would be useful to quantify the degree of rup-
ture asymmetry in various cases by measuring (as done in Fig. 2.7 and related phase dia-
grams) the potency values that are generated in the positive and negative directions.
Wrinkle-like ruptures have remarkable dynamic properties that may be relevant to
many issues of earthquake and fault mechanics [e.g., Ben-Zion, 2001]. One important
unresolved issue is the range of realistic conditions (if any) for which wrinkle-like rup-
tures can produce opening modes and the associated dynamic process. Clarifying rupture
properties over the range of parameters for which our calculations develop large numeri-
cal oscillations (white area in Fig. 2.12) will provide a better understanding of the transi-
tion between crack-like and pulse-like modes of rupture. Finally, it is important to per-
form additional in situ tests of signals predicted for wrinkle-like ruptures as done in the
seismological and geological studies of Rubin and Gillard [2000], McGuire et al. [2002],
Lewis et al. [2005] and Dor et al. [2006a,b].
36
Chapter 3
Energy Partition During Dynamic Rupture
Process
3.1 Introduction
The occurrence of dynamic rupture converts elastic strain energy that is stored in the
medium to other forms of energy. These consist primarily of kinetic energy that is radi-
ated away from the rupture, dissipative energy near the failure surface in the form of heat,
and dissipation in the bulk due to the creation of new surface area and inelastic strain
[e.g., Broberg, 1999; Aki and Richards, 2002; Ben-Zion, 2003]. However, the fractional
values of the components and other details associated with the energy partition process
are not well understood. The properties of dynamic ruptures and energy partition have
fundamental implications for many issues of both basic and applied science. On the basic
side, a better understanding of the energy partition process can provide guidance for the
development of an improved quantitative theory for dynamic rupture. Important applica-
tions include constraints for seismic hazard associated with the energy radiated during
earthquake ruptures and expected frictional heat in the vicinity of earthquake faults. The
overall goal of this study is to obtain a quantitative understanding of the relations be-
tween properties of dynamic ruptures and the energy partition process under different
37
conditions. As an initial step in that direction, we focus in the present study on a two-
dimensional plane-strain homogeneous framework with a planar frictional interface gov-
erned by a rate- and state-dependent (RSD) friction law.
Numerous works have shown that rupture on a frictional interface in a homogeneous
solid can occur in either a crack-like mode or a pulse-like mode [e.g., Das and Kostrov,
1988; Beeler and Tullis, 1996; Beroza and Mikumo, 1996; Zheng and Rice, 1998; Niel-
sen et al., 2000]. In the former case the slipping region expands continuously until rup-
ture terminates, whereas in the latter case only a small portion of the interface slips at any
given time. The crack mode is favored by relatively high and smooth initial background
stress, relatively weak velocity dependence of friction, and little or no coupling of slip to
changes of normal stress, while the pulse mode is favored by the opposite set of condi-
tions [e.g., Ben-Zion, 2001; Rice, 2001]. Recent numerical simulations have shown that
for ranges of conditions there is a third transitional mode of rupture consisting of a train
of pulses [Coker et al., 2005; Lapusta, 2005; Shi and Ben-Zion, 2006].
Another topic of considerable interest is the speed of rupture propagation [e.g., Rosa-
kis et al., 1999; Rosakis, 2002]. It is generally understood that increasing the initial shear
stress on a frictional interface toward the yield strength leads to a transition from sub-
shear to supershear propagation [Burridge, 1973; Andrews, 1976; Das and Aki, 1977;
Day, 1982; Xia et al., 2004]. Since increasing initial shear stress favors the crack-like
rupture mode, one would expect supershear propagation to be associated primarily with
crack-like rather than pulse-like mode of rupture. Inversions of seismic data suggest that
portions of earthquake ruptures are associated occasionally with supershear velocities
38
[e.g., Archuleta, 1984; Bouchon and Vallée, 2003]. However, it is generally accepted that
earthquake rupture speeds are typically about 75% of the shear wave velocity [e.g., Ben-
Zion, 2003] and that earthquakes tend to propagate as pulses [e.g., Heaton, 1990].
Rupture on a frictional interface belongs to the class of non-linear problems for which
there is strong coupling between small-scale processes and the overall large-scale proper-
ties of the response. Rice [1993], Ben-Zion and Rice [1995, 1997], Rice et al. [2001] and
others emphasized the coupling between small-scale features of the friction and large-
scale aspects of the slip process. Recent studies have also pointed out that the mode,
speed and other properties of simulated ruptures can be strongly affected by the em-
ployed nucleation procedure [Shi and Ben-Zion, 2006; Festa and Vilotte, 2006; Ben-
Zion, 2006].
In the current chapter we examine properties of dynamic rupture and energy partition
by performing two-dimensional finite element simulations of rupture along a planar fric-
tional interface between identical isotropic elastic materials. The obtained mode (crack,
pulse, train of pulses), speed (subshear or subshear) and energy partition associated with
the rupture process are shown to depend on the assumed friction and nucleation parame-
ters and the initial shear stress. In Section 3.2.2 we describe the model setup and the RSD
constitutive law employed in the simulations. In Section 3.2.3 we present results associ-
ated with four cases having different rupture properties and partition of energy between
the radiation and heat on the interface. The results are discussed and summarized in Sec-
tion 3.2.4.
39
3.2 Model Setup
3.2.1 Finite Element Formulation
A two-dimensional in-plane dynamic rupture model of two blocks separated by a pla-
nar interface S
int
is studied numerically using the finite element method. A Lagrangian
description of small-strain deformation is employed and the deformed state with initial
stresses is taken as the reference configuration. The principle of virtual work is written as
(Needleman, 1987)
€
σ :δεdA
A
∫
− T⋅δΔudS
S
int
∫
= T⋅δudS
S
ext
∫
− ρ
∂
2
u
∂t
2
δudA
A
∫
, (3.1)
where t is time, σ the stress tensor, ε is the strain tensor, u is the displacement vector, T is
the traction vector, Δu is the displacement jump across the interface S
int
and σ:δε denotes
σ
ij
δε
ij
. The model area, internal interface between the two blocks and external boundaries
in the reference configuration are, respectively, A, S
int
and S
ext
. With subscripts ‘n’ and ‘s’
denoting in-plane directions normal and tangential to the interface S
int
, respectively, T
n
and T
s
are the normal and tangential components of the traction vector T along S
int
.
The finite element discretization is based on linear displacement triangular elements
that are arranged in a ‘cross-triangle’ quadrilateral pattern. When the discretization is
substituted into the principle of virtual work of Eq. (3.1), and the integrations are carried
out, the discretized equations of motion become
€
M
∂
2
U
∂t
2
=R, (3.2)
where U is the vector of nodal displacement, M is the mass matrix and R is the nodal
force vector consisting of contributions from the area elements and the cohesive surfaces.
40
To achieve higher accuracy and computational efficiency with the employed explicit time
integration procedure, a lumped mass matrix is used in Eq. (3.2) instead of the consistent
mass matrix [Krieg and Key, 1973]. Four-point Gaussian integration is used along the
interface S
int
. The time integration of the governing equation (3.2) is based on a central-
difference scheme which is a member of the Newmark β-method [Belytschko et al.,
1976].
3.2.2 Model Configurations
Figure 3.1: (a) Geometry of the two-dimensional plain-strain model and the initial load-
ing conditions. The dashed bold line delineates the external boundaries. The frictional
interface S
int
is along x
2
= 0. The nucleation zone S
nucl
is a confined central portion of S
int
denoted by a thick solid line. (b) The ‘cross-triangle’ quadrilateral mesh pattern em-
ployed to discretize the model. See text for additional details.
The geometry of the two-dimensional in-plane model and its initial loading condi-
tions are illustrated in Fig. 3.1(a). A Cartesian coordinate system is used with the x
1
− x
2
plane being the plane of deformation. The interface S
int
is along x
2
= 0. The model has a
(a) (b)
Τ
0
−Σ
0
Τ
0
−Σ
0
x
2
nucleation zone
S
nucl S
int
S
ext
x
1
fine uniform
coarse uniform gradual
41
square geometry of −5 m ≤ x
1
≤ 5 m and −5 m ≤ x
2
≤ 5 m with the origin of the coordi-
nate system (x
1
, x
2
) = (0, 0) being the center. Uniform compressive normal stresses
€
σ
11
0
=σ
22
0
=−Σ
0
and shear stress
€
σ
12
0
=T
0
are applied at the external boundaries S
ext
. The
region S
nucl
is a confined segment of S
int
, symmetric about the origin, where we impose
nucleation process that initiates dynamic rupture events (see Sections 3.2.4 and 3.2.5).
Fig. 3.1(b) illustrates the triangulated-quadrilateral mesh pattern employed in this
study. The mesh consists of a central region with fine discretization in the form of uni-
form rectangles each of which is further divided into four triangles. The central fine mesh
is extended outwards first by a gradual mesh and then by a coarse uniform mesh near the
external boundaries.
The fine-mesh region in the model has a size of 4 m × 4 m (−2 m ≤ x
1
≤ 2 m, −2 m ≤
x
2
≤ 2 m) and is discretized into 160,000 quadrilaterals (640,000 triangular elements).
Each quadrilateral is a square with a size of 0.01 m × 0.01 m. The magnitude of the initial
normal stress Σ
0
is set to be 100 MPa for all the simulations in this chapter.
The materials across the interface S
int
are identical isotropic elastic media with
Young’s modulus E, Poisson’s ratio ν and density ρ. For plain-strain problems, the speeds
of the P wave (α), S wave (β) and Rayleigh wave (c
R
) are
€
α =
1−ν ( )E
ρ 1+ν ( ) 1−2ν ( )
, (3.3)
€
β =
E
2ρ 1+ν ( )
, (3.4)
€
c
R
=2β 1−γ
2
2−
4
3
γ + R + M
3
+ R− M
3
−1 2
, (3.5)
42
where
€
R =2 27−90γ +99γ
2
−32γ
3
( )
27,
€
M =4 1−γ
( )
2
11−62γ +107γ
2
−64γ
3
( )
27 and
€
γ = β α ( )
2
.
The elasticity parameters of the model are chosen to represent Homalite-100: E = 5.3
GPa, ν = 0.35 and ρ = 1246 kg/m
3
. With these values the elastic wave speeds are α =
2613 m/s, β = 1255 m/s and c
R
= 1174 m/s. The size of the time step in our simulations,
determined on the fly based on the sliding condition, ranges from 10
−6
× min(ds) to 10
−2
×
min(ds) , where min(ds) is the shortest edge of all the triangular elements. The required
very small fractions of the Courant time step are typical to simulations of dynamic rup-
tures on an interface governed by RSD friction [e.g., Ben-Zion and Rice, 1997; Lapusta
et al., 2000].
We examine the rupture dynamics only in the central region of the model discretized
by the fine uniform mesh. The outer part of the model is used to accommodate the stress
waves that are radiated from the nucleation zone and the ensuing dynamic rupture. Each
calculation is terminated before the rupture propagation reaches the edges of the fine-
mesh region or alteration of the stress field reaches the external boundaries S
ext
. Therefore
there is zero energy input from the external boundaries S
ext
during the entire simulation
process, i.e., from the start of the nucleation process to the termination of the calculation.
As a consequence, the first term on the right-hand side of Eq. (3.1) vanishes during our
simulations.
43
3.2.3 Constitutive relations
The current model employs two constitutive relations: Hookean elasticity relating
stress and strain in the volume and a general form of RSD friction law characterizing the
frictional response of the interface S
int
during sliding. Following Coker et al. [2005], we
adopt the following expression of RSD friction coefficient µ proposed by Povirk and
Needleman [1993]
€
µθ
0
,Δ˙ u
slip
( )
= gθ
0
( )
1+Δ˙ u
slip
V
0
( )
1 m
(3.6)
with
€
gθ
0
( )
=
µ
d
+ µ
s
− µ
d
( )
exp − L
0
θ
0
( )
V
1
( )
p
[ ]
1+ L
0
θ
0
( )
V
0
( )
1 m
, (3.7)
where θ
0
is a state variable with evolution given by
€
˙
θ
0
= B 1− θ
0
Δ˙ u
slip
( )
L
0 [ ]
. (3.8)
In Eqs (3.6)-(3.8), µ
s
and µ
d
are the nominal static and dynamic coefficients of fric-
tion, L
0
is the characteristic slip distance over which the friction coefficient evolves to a
steady-state value, V
0
and V
1
are reference slip velocities and B is a constant.
The steady-state value µ
ss
of the friction coefficient, obtained by setting
€
˙
θ
0
= 0 (for
which
€
θ
0
= L
0
Δ˙ u
slip
), is
€
µ
ss
= µ
d
+ µ
s
− µ
d
( )exp −Δ˙ u
slip
V
1 ( )
p
[ ]
. (3.9)
The initial value of the state variable θ
0
is set to be
€
θ
0
0
( )
=hL
0
V
1
, (3.10)
where h is a constant.
44
Setting
€
Δ˙ u
slip
= 0 in Eq. (3.6) and replacing θ
0
with Eq. (3.10), we obtain an expres-
sion for the apparent static coefficient of friction µ
as
as
€
µ
as
=
µ
d
+ µ
s
− µ
d
( )
exp−h
−p
[ ]
1+hV
1
V
0
( )
1 m
. (3.11)
Unless otherwise specified, the values of frictional parameters that are used in our
simulations are µ
s
= 0.6, µ
d
= 0.5, L
0
= 0.2 mm, B = 4.6, m = 5 and p = 1.2. Some of these
values were found in previous studies to characterize experimental observations of
Homalite sliding on Homalite [Rosakis et al., 1999, 2000; Rosakis, 2002; Samudrala et
al., 2002; Coker et al., 2005].
The frictional response associated with Eqs (3.6)-(3.8) is illustrated in Fig. 3.2. An
abrupt change of slip rate produces an instantaneous effect with the same sign as the
change of slip rate, followed by an evolutionary effect over a characteristic slip distance
with an opposite sign. This behavior was shown to characterize the frictional response of
different types of rocks and other materials [Dieterich, 1979, 1981; Ruina, 1983;
Dieterich and Kilgore, 1996]. In the employed formulation, the parameter V
0
controls the
amplitude of direct effect of the response of the apparent friction coefficient µ
app
, defined
as T
s
/ T
n
, to changes of slip rate. Increasing V
0
reduces the direct effect (Fig. 3.2a). The
parameter V
1
is a threshold slip rate for which µ
app
approaches the nominal dynamic level
µ
d
. Increasing V
1
spreads the rate sensitivity of the frictional response over a broader
range of slip rates (Fig. 3.2b). The steady-state coefficient of friction µ
ss
as a function of
slip rate is shown for different values of V
1
in Fig. 3.2(c). In the simulations below we use
V
0
= 2000 m/s in conjunction with V
1
= 0.1 m/s or 0.65 m/s. With these parameters, the
45
response of µ
app
to abrupt changes of slip rate is small (Fig. 3.2d) and the employed RSD
friction mimics the behavior of slip-weakening friction with a characteristic slip distance
L
0
when the slip rate is sufficiently large (
€
Δ˙ u
slip
>> V
1
). To estimate the appropriate nu-
merical grid size and critical nucleation length, we use below results based on a linear
stability analysis for a classical slip-weakening model with peak strength µ
as
σ
n
, dynamic
strength µ
d
σ
n
and characteristic slip distance L
0
.
Figure 3.2: Effects of parameters of the employed rate- and state-dependent friction. (a)
The parameter V
0
controls the amplitude of the direct response of the apparent coefficient
of friction µ
app
to abrupt changes of slip rate. Increasing V
0
reduces the direct effect. (b)
The parameter V
1
serves as a threshold slip rate for which µ
app
approaches the nominal
dynamic level µ
d
. Increasing V
1
increases the range of slip rates for which the frictional
response is sensitive. (c) The steady-state coefficient of friction as a function of slip rate
for different values of V
1
. (d) Evolution of µ
app
in response to abrupt changes of slip rate
under constant compressive normal stress for two combinations of friction parameters
used in the simulations.
0 5 10 15 20 25 30 35 40
0.5
0.52
0.54
0.56
0.58
0.6
slip rate (m/s)
µ
ss
V
= 26 m/s
V
= 15 m/s
V
= 5 m/s
V
= 0.65 m/s
V
= 0.1 m/s
0 2 4 6 8 10 12
0.45
0.5
0.55
0.6
0.65
Δu
slip
/L
0
µ
app
V
= 0.65 m/s
V
= 0.1 m/s
1 m/s 30 m/s
(c)
(d)
0 2 4 6 8 10
0.4
0.45
0.5
0.55
0.6
0.65
Δu
slip
/L
0
µ
app
V
= 100 m/s
V
= 300 m/s
V
= 2000 m/s
(a)
0 2 4 6 8 10
0.4
0.45
0.5
0.55
0.6
0.65
Δu
slip
/L
0
µ
app
V
= 26 m/s
V
= 10 m/s
V
= 0.1 m/s
(b)
100 m/s
10 m/s 30 m/s 100 m/s 10 m/s 30 m/s 100 m/s
V
= 26 m/s V
= 100 m/s
V
= 2000 m/s
46
For sliding along an interface governed by a simple linear slip-weakening friction, the
spatial region behind the rupture tip where the strength degradation occurs is given ap-
proximately [Rice, 1980; Ben-Zion, 2003] by
€
X
c
=cGL
0
σ
n
µ
as
− µ
d
( ) [ ]
, (3.12)
where c is a dimensionless constant of about 1−3. Using normal stress σ
n
= Σ
0
= 100
MPa, shear modulus G = 1.963 GPa, µ
as
− µ
d
≈ 0.1 and L
0
= 0.2 mm we have X
c
≈ 0.1 m.
The employed grid spacing of 0.01 m in the fine mesh region is 10 times smaller than this
estimate.
Frictional sliding experiments with variable normal stress show [Prakash and Clifton,
1993; Prakash 1998] that the shear strength responds gradually to abrupt changes of nor-
mal stress. In the simulations performed in this study the normal traction along the inter-
face S
int
remains constant, so the Prakash-Clifton type response, incorporated in the fric-
tion formulation of Coker et al. [2005], is not included. Simulations with model realiza-
tions incorporating material contrast and faults with geometrical irregularities require
consideration of the Prakash-Clifton response.
In addition to the frictional behavior of the interface discussed above, the calculations
incorporate a linear elastic stiffness between material points directly opposing each other
across S
int
characterized by two stiffness parameters C
n
and C
s
[Povirk and Needleman,
1993]. Depending on the amount of the normal displacement jump across S
int
, the follow-
ing two situations can exist. (I) If Δu
n
> Σ
0
/C
n
, the interface is in tension and tractions
along the associated region vanish. (II) If Δu
n
≤ Σ
0
/C
n
, the two sides of the interface are in
contact and the frictional law is invoked. In the latter case,
47
€
˙
T
n
=−C
n
Δ˙ u
n
, (3.13)
€
˙
T
s
=−C
s
Δ˙ u
s
− sgn T
s
( )Δ˙ u
slip [ ]
, (3.14)
where
€
Δ˙ u
n
and
€
Δ˙ u
s
are, respectively, velocity differences of velocity in the normal and
tangential directions between material points that are on opposite sides of the interface in
the initial configuration. Note that
€
Δ˙ u
slip
is the magnitude of the pure frictional slip rate
specified by the RSD friction law (
€
Δ˙ u
slip
≥ 0). The jump of the tangential particle velocity
across the interface
€
Δ˙ u
s
is the sum of two parts
€
Δ˙ u
s
=Δ˙ u
s
el
+Δ˙ u
slip
, (3.15)
where
€
Δ˙ u
s
el
=−
˙
T
s
C
s
. In the simulations performed in this study
€
Δ˙ u
n
= 0, so only situa-
tion (II) is applicable here.
Coker et al. [2005] examined frictional sliding modes along an interface between two
identical Homalite-100 plates subjected to shear impact loading. Based on comparisons
between numerical simulations and laboratory experiments they concluded that C
n
= 300
GPa/m and C
s
= 100 GPa/m provide a reasonable agreement with the experimental data.
These values will also be used in all our simulations in this chapter. An important pa-
rameter that governs the rupture speed with slip-weakening friction is the strength excess
parameter S [e.g., Andrews, 1976; Das and Aki, 1977] defined as
€
S = µ
as
σ
n
−T
0
( )
T
0
− µ
d
σ
n
( )
, (3.16)
where µ
as
σ
n
is the peak strength, T
0
is the initial shear stress, and µ
d
σ
n
is the dynamic
strength level of the interface.
48
3.2.4 Nucleation Procedure
In addition to the friction law and stress heterogeneity, the type of employed nuclea-
tion procedure and related parameters can also strongly affect the rupture mode. Ruptures
in our model are nucleated by gradually increasing uniformly over time t
nucl
the shear
stress in a confined central region S
nucl
with width L
nucl
(Fig. 3.1a). The final shear stress
level at time t
nucl
is typically about 2 MPa above the strength of the interface with t
nucl
=
20 µs. The RSD friction law is temporarily deactivated in the nucleation zone S
nucl
during
the stage of rupture nucleation (0 ≤ t ≤ t
nucl
).
Based on the quasi-static analysis on the instability of expanding slip patch, Uenishi
and Rice [2003] proposed a critical length for linear slip-weakening model that depends
only on the shear modulus and the slip-weakening rate:
€
L
c
=1.158G
*
γ , (3.17)
where G
*
= G/(1−ν) for mode II and γ = (µ
as
− µ
d
)σ
n
/L
0
is the weakening rate. With our
employed parameters Eq. (3.17) gives L
c
= 0.074 m, a value smaller than the 0.1 m esti-
mate of degradation/nucleation size associated with Eq. (3.12). In the simulations of Sec-
tion 3.2.3, we choose L
nucl
to be 0.2 m for most cases and a larger value of 0.4 m for a
case that produces a train of pulses.
3.2.5 Energy Calculations
The entire simulation process can be separated into two stages: the first stage is rup-
ture nucleation and the second stage is spontaneous rupture propagation. Since we termi-
nate the calculations before the stress waves reach the external boundaries, there is zero
49
energy input from the tractions at S
ext
. Then from the principle of virtue work, the power
balance equation can be written as
€
σ
ij
˙ ε
ij
dA
A
∫
− T
i
⋅Δ˙ u
i
dS
S
int
−S
nucl
∫
− T
i
⋅Δ˙ u
i
dS
S
nucl
∫
+ ρ˙ ˙ u
i
˙ u
i
dA
A
∫
= 0. (3.18)
Note that all the tensors and vectors are written in a component form and the contri-
bution from the nucleation zone S
nucl
is nonzero only during the first stage (i.e., rupture
nucleation) and vanishes for t > t
nucl
.
Integration of each term in Eq. (3.18) over time gives following energy components:
(i) Change of elastic strain energy in the volume:
€
W
vol
el
= σ
ij
˙ ε
ij
dA
A
∫ [ ]
dt
0
t
∫
. (3.19)
(ii) Energy input by the nucleation procedure:
€
W
app
= T
i
⋅Δ˙ u
i
dS
S
nucl
∫
[ ]
dt
0
t
∫
. (3.20)
(iii) Kinetic energy:
€
W
kin
=
1
2
ρ˙ u
i
˙ u
i
dA
A
∫
. (3.21)
(iv) Change of interface energy:
€
W
int
=− T
i
⋅Δ˙ u
i
dS
S
int
−S
nucl
∫
[ ]
dt
0
t
∫
. (3.22)
The interface energy (3.22) can further be written as
€
W
int
=− T
i
⋅Δ˙ u
i
dS
S
int
−S
nucl
∫
[ ]
dt
0
t
∫
€
=− T
s
⋅Δ˙ u
s
dS
S
int
−S
nucl
∫
[ ]
dt
0
t
∫
− T
n
⋅Δ˙ u
n
dS
S
int
−S
nucl
∫
[ ]
dt
0
t
∫
. (3.23)
50
Substitution of Eq. (3.15) into the first term of Eq. (3.23) gives
€
− T
s
⋅Δ˙ u
s
dS
S
int
−S
nucl
∫
[ ]
dt
0
t
∫
€
=− T
s
⋅Δ˙ u
s
el
dS
S
int
−S
nucl
∫
[ ]
dt
0
t
∫
− T
s
⋅Δ˙ u
slip
dS
S
int
−S
nucl
∫
[ ]
dt
0
t
∫
€
= T
s
⋅
˙
T
s
C
s
dS
S
int
−S
nucl
∫
[ ]
dt
0
t
∫
− T
s
⋅Δ˙ u
slip
dS
S
int
−S
nucl
∫
[ ]
dt
0
t
∫
. (3.24)
The second term in Eq. (3.24) corresponds to dissipative frictional heat
€
W
heat
=− T
i
⋅Δ˙ u
slip
dS
S
int
−S
nucl
∫
[ ]
dt
0
t
∫
, (3.25)
while the elastic part of the interface energy is
€
W
int
el
= T
s
⋅
˙
T
s
C
s
dS
S
int
−S
nucl
∫
[ ]
dt
0
t
∫
− T
n
⋅Δ˙ u
n
dS
S
int
−S
nucl
∫
[ ]
dt
0
t
∫
. (3.26)
Using the above definitions, the balance equation for energy changes can be written
as
€
W
vol
el
+W
int
+W
kin
=W
app
, (3.27)
or
€
W
el
+W
heat
+W
kin
=W
app
, (3.28)
where
€
W
el
=W
vol
el
+W
int
el
is the total change of elastic energy.
3.3 Simulation results
By varying the initial shear stress T
0
which is related to the strength excess parameter
S, the velocity dependence of the friction law (friction parameter V
1
) and the width of the
nucleation zone L
nucl
, we obtain the following four cases with different rupture modes:
Case I – supershear crack-like rupture, Case II – subshear crack-like rupture, Case III –
51
subshear single rupture pulse, and Case IV – supershear train of pulses. The parameters
used to obtain these four cases are:
Case I: T
0
= 58 MPa, S = 0.17, V
1
= 0.1 m/s and L
nucl
= 0.2 m;
Case II: T
0
= 54 MPa, S = 1.35, V
1
= 0.1 m/s and L
nucl
= 0.2 m;
Case III: T
0
= 54 MPa, S = 1.35, V
1
= 0.65 m/s and L
nucl
= 0.2 m;
Case IV: T
0
= 54 MPa, S = 1.35, V
1
= 0.65 m/s and L
nucl
= 0.4 m.
Fig. 3.3 shows time histories of the shear stress and normalized slip, along with the
functional dependence of the shear stress on normalized slip at a fixed point x
1
= 1.0 m on
the interface S
int
for all four cases. The curves of stress vs. normalized slip clearly exhibit
a slip-weakening-like behavior that is expected based on the employed set of frictional
parameters. In each case, the rupture propagates all the way from the edge of the nuclea-
tion zone to the edge of the fine-mesh region. The amount of slip generated at that dis-
tance is largest for the supershear crack (Fig. 3.3a) and smallest for the single pulse (Fig.
3.3c). The subshear crack (Fig. 3.3b) and the train of pulses (Fig. 3.3d) have intermediate
amounts of slip. The corresponding amount of slip in Fig. 3.3(b) is about 20% less than in
Fig. 3.3(a) and for the train of pulses (Fig. 3.3d) it is about half of that in Fig. 3.3(b). The
amount of slip for the single pulse (Fig. 3.3c) is nearly an order of magnitude smaller
than for the supershear crack. The small steps in the curve of slip vs. time of Fig. 3.3(d)
are footprints of the train of pulses. We note that for case III with a rupture pulse (Fig.
3.3c), the shear stress recovers nearly to the initial level after the passage of the pulse. A
similar situation holds for case IV with the train of pulses (Fig. 3.3d), but after a larger
amount of slip than is shown.
52
Figure 3.3: Time histories of slip and shear stress (left columns) and functional depend-
ence of the shear stress on slip (right columns) at the point x
1
= 1.0 m on the interface for
(a) Case I with supershear crack-like rupture, (b) Case II with subshear crack-like rup-
ture, (c) Case III with subshear single pulse, and (d) Case IV with supershear train of
pulses.
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8
50
52
54
56
58
60
Time (ms)
Shear stress (MPa)
0
5
10
15
20
25
Δu
slip
/L
0
0 0.5 1 1.5 2 2.5 3
50
52
54
56
58
60
Δu
slip
/L
0
Shear stress (MPa)
(a)
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6
50
52
54
56
58
60
Time (ms)
Shear stress (MPa)
0
4
8
12
16
20
Δu
slip
/L
0
0 0.5 1 1.5 2 2.5 3
50
52
54
56
58
60
Δu
slip
/L
0
Shear stress (MPa)
(b)
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6
50
52
54
56
58
60
Time (ms)
Shear stress (MPa)
0
0.5
1
1.5
2
2.5
3
Δu
slip
/L
0
0 0.5 1 1.5 2 2.5 3
50
52
54
56
58
60
Δu
slip
/L
0
Shear stress (MPa)
(c)
0 0.2 0.4 0.6 0.8 1 1.2
50
52
54
56
58
60
Time (ms)
Shear stress (MPa)
0
2
4
6
8
10
Δu
slip
/L
0
0 0.5 1 1.5 2 2.5 3
50
52
54
56
58
60
Δu
slip
/L
0
Shear stress (MPa)
(d)
53
Figure 3.4: Rupture front speeds V
R
in the positive x
1
direction versus time for Cases I-
IV. The speeds of the P wave (α), S wave (β) and
€
2β are indicated by the dashed lines.
Fig. 3.4 shows the values of the rupture front speed V
R
in the +x
1
direction for all four
cases. The positions of the rupture fronts in both directions (+x
1
and − x
1
) are recorded at
various times and V
R
is then calculated by a seven-point progressive least-squares fit on
these position data. Since we use a model with the same elastic properties across the in-
terface S
int
, the rupture characteristics in the opposite directions are identical. In Case I,
the rupture front quickly accelerates beyond
€
2β after leaving the nucleation region and
then gradually approaches the P wave speed α. In Case II, V
R
follows closely the
Rayleigh wave speed c
R
. The self-healing pulse in Case III travels at a speed of about
0.7c
R
after some initial transients. In Case IV, the rupture propagates initially after the
nucleation in a crack-like mode with V
R
between β and
€
2β . The middle portion of the
sliding crack heals quickly producing two separate slip patches that propagate in the op-
posite directions with accelerating rupture front speeds. Each slip patch then splits gradu-
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
0
0.5
1.0
1.5
2.0
2.5
3.0
Time (ms)
Rupture front speed V
R
(km s
-1
)
β
β
α
c
R
Case I
Case II
Case III
Case IV
2
54
ally into a train of pulses and additional pulses are generated near the end of the existing
train of pulses adding to the overall length of the train. During this process the rupture
front speed approaches a value that is midway between
€
2β and α. These and other de-
tails of the rupture behavior are seen more clearly in subsequent figures. The oscillations
of the rupture front speeds in Fig. 3.4 decrease when we reduce the mesh spacing. How-
ever, the estimated average values of V
R
and other features discussed in this work remain
essentially the same.
The profiles of frictional slip rate and shear traction along the interface at time 0.7 ms
for the four cases are shown in Fig. 3.5. The slip rates (red) exhibit high-frequency oscil-
lations that are especially pronounced in Cases I and IV having supershear propagation
speeds. As shown in the insets, the oscillations are well resolved numerically. Each wig-
gle in the oscillatory profile consists of at least 14 points for Case I (Fig. 3.5a) and at least
8 points in Case II (Fig. 3.5b). The pulse in Case III (Fig. 3.5c) has more than 24 points
consistently and each single pulse in Case IV (Fig. 3.5d) has about 12 points. The shear
stress profiles (blue) show high-frequency oscillations behind the expanding sliding areas
in Cases III and IV that are associated with pulses. The high-frequency oscillations of slip
rates and shear stress may be produced by the relatively rapid nucleation process, the
presence of the elastic interface stiffness parameter C
s
, the large dynamic shear stress
drop and the large reduction of the friction coefficient at relatively low sliding rate. We
repeated the calculations with both finer grid (0.005 m spacing in the fine-mesh region)
and coarser grid (0.02 m spacing in the fine-mesh region) for all four cases and found that
the general features of the rupture properties and magnitudes of all the energy compo-
55
nents remain the same, although perfect wiggle-to-wiggle match is not achieved with dif-
ferent mesh resolutions.
Figure 3.5: Profiles of frictional slip rate (red) and shear traction (blue) along the inter-
face at time 0.7 ms for (a) Case I, (b) Case II, (c) Case III, and (d) Case IV. The insets
show close-up views of the wiggles in the slip-rate profiles.
In the crack-like Cases I and II, the slip rate is positive everywhere behind the propa-
gating rupture fronts (Figs 3.5a-b). In Case III with a pulse, the slip rate is nonzero only
in narrow regions behind the fronts (Fig. 3.5c). In Case IV with a train of pulses, the slip
rate alternates between positive and zero values in a broadening region behind the fronts
(Fig. 3.5d). The shear stress profiles for Cases I and II show the smooth characteristics of
crack-like ruptures, while the shear stress profile for Case III with a pulse has an addi-
-2 -1.5 -1 -0.5 0 0.5 1 1.5 2
50
52
54
56
58
60
Shear stress (MPa)
0
4
8
12
16
20
Slip rate (m/s)
-2 -1.5 -1 -0.5 0 0.5 1 1.5 2
50
52
54
56
58
60
Position along interface (m)
Shear stress (MPa)
0
20
40
60
80
100
Slip rate (m/s)
-2 -1.5 -1 -0.5 0 0.5 1 1.5 2
50
52
54
56
58
60
Position along interface (m)
Shear stress (MPa)
0
10
20
30
40
50
60
Slip rate (m/s)
-2 -1.5 -1 -0.5 0 0.5 1 1.5 2
50
52
54
56
58
60
Position along interface (m)
Shear stress (MPa)
0
4
8
12
16
20
Slip rate (m/s)
(a) Case I (b) Case II
(c) Case III
(d) Case IV
Position along interface (m)
56
tional localized structure. The stress profiles for Cases III and IV have persisting large-
scale oscillations in the nucleation zones with amplitudes that decay with time.
Figure 3.6: Time histories of various energy components for (a) Case I, (b) Case II, (c)
Case III and (d) Case IV. The insets illustrate the approximate partition of the total
change of elastic energy between the kinetic energy W
kin
and heat energy W
heat
.
Fig. 3.6 gives the evolution of various energy components for each case with insets
showing the partition of the total change of elastic energy between the kinetic energy (ra-
diation) W
kin
and heat energy (dissipation) W
heat
. For the same rupture distance, Case I
with a supershear crack has the largest amount of kinetic energy, while Case III with a
single subshear pulse has the smallest amount of radiated energy in terms of both the
magnitude and percentage value. As shown in Fig. 3.7, the mismatch in percentage be-
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8
−1
−0.8
−0.6
−0.4
−0.2
0
0.2
0.4
0.6
0.8
1
Time (ms)
0 0.6 1.0 1.4 2.0
−0.8
−0.6
−0.4
−0.2
0
0.2
0.4
0.6
0.8
Time (ms)
0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6
−0.08
−0.06
−0.04
−0.02
0
0.02
0.04
0.06
0.08
Time (ms)
0 0.2 0.4 0.6 0.8 1.0
−0.25
−0.2
−0.15
−0.1
−0.05
0
0.05
0.1
0.15
0.2
0.25
Time (ms)
Energy (x10
6
J/m)
Energy (x10
6
J/m)
Energy (x10
6
J/m)
Energy (x10
6
J/m)
0.2 0.4 0.8 1.2 1.6 1.8
W
vol
el
Wapp
W
heat
W
int
el
W
kin
(94%)
(6%)
W
heat
W
kin
W
vol
el
Wapp
W
heat
W
int
el
W
kin
W
vol
el
Wapp
W
heat
W
int
el
W
kin
W
vol
el
Wapp
W
heat
W
int
el
W
kin
(97%)
(3%)
W
heat
W
kin
(99%)
(1%)
W
heat
W
kin
(98%)
(2%)
W
heat
W
kin
(a) Case I (b) Case II
(c) Case III
(d) Case IV
57
tween the sum of energy components
€
W
vol
el
+W
int
el
+W
heat
+W
kin
and
€
W
app
is generally
less than 0.4% except for Case I where it is ~0.6%. Therefore the energy balance equa-
tion (3.28) is well satisfied numerically.
Figure 3.7: The mismatch in percentage between the energies
€
W
vol
el
+W
int
el
+W
heat
+W
kin
and
€
W
app
for Cases I-IV.
Figs 3.8-3.11 show the spatial distributions of the shear stress σ
12
and snapshots of
slip profiles for all four cases. The characteristic Mach cones of supershear rupture are
seen clearly in the stress contours of Cases I (Fig. 3.8a) and IV (Fig. 3.11a). The gradual
distributions of stress behind the rupture fronts in Figs 3.8(a) and 3.9(a) imply that slip
occurs continuously along the interface behind the propagating rupture fronts. This is
confirmed by the evolutions of the slip profiles along the interface (Figs 3.8b and 3.9b)
that follow crack-like circular shapes. In contrast, the localized concentrations of stress
0 0.5 1 1.5 2
0
0.2
0.4
0.6
0.8
1.0
1.2
1.4
1.6
1.8
2
Time (ms)
Energy Mismatch in Percentage
Case I
Case II
Case III
Case IV
58
behind the rupture fronts in Figs 3.10(a) and 3.11(a), and the narrow propagating slip
fronts with overall triangular profiles in Figs 3.10(b) and 3.11(b) are indicative features
of pulses. Manighetti et al. [2005] analyzed slip profiles of several tens of large earth-
quakes and concluded that most of them are characterized by triangular shapes.
For the four cases presented above, the nominal dynamic coefficient µ
d
is fixed at 0.5,
a large portion of the nominal static coefficient 0.6. Recent experiments with rock friction
suggest that the reduction of the friction coefficient can be significant at high seismic slip
rates [e.g., Tsutsumi and Shimamoto, 1997; Goldsby and Tullis, 2002; Di Toro et al.,
2004]. To illustrate some tendencies associated with lower dynamic friction we show in
Fig. 3.12 results based on µ
d
= 0.3 and all other parameters same as in Case II. With the
strength excess parameter S = 0.17, the rupture propagates in a supershear crack-like
fashion (Fig. 3.12c) similar to Case I, but with larger slip rate (Fig. 3.12a) and higher rup-
ture front speed (Fig. 3.12b). The fraction of the radiated energy (Fig. 3.12d) in this case
is 12%, which is considerably higher than all four previous cases. However, we note that
the absolute magnitude of the generated heat energy is also higher than all the previous
cases, since the increase in the slip rate has a larger effect than the decrease in the sliding
friction level.
59
Figure 3.8: (a) The spatial distribution of shear stress σ
12
at time 0.5 ms and (b) snapshots
of the slip profiles along the interface at various times for Case I with supershear crack-
like rupture.
x
1
(m)
x
2
(m)
-2 -1 0 1 2
-2
-1
0
1
2
56
58
60
62
64
66
68
σ
12
-2 -1.5 -1 -0.5 0 0.5 1 1.5 2
0
1
2
3
4
5
6
Position along interface (m)
Slip (mm)
0.15 ms
0.3 ms
0.45 ms
0.6 ms
0.75 ms
(a)
(b)
60
Figure 3.9: (a) The spatial distribution of shear stress σ
12
at time 0.7 ms and (b) snapshots
of the slip profiles along the interface at various times for Case II with subshear crack-
like rupture.
σ
12
(a)
(b)
x
1
(m)
x
2
(m)
-2 -1 0 1 2
-2
-1
0
1
2
52
54
56
58
σ
12
-2 -1.5 -1 -0.5 0 0.5 1 1.5 2
0
0.5
1
1.5
2
2.5
3
3.5
4
Position along interface (m)
Slip (mm)
0.3 ms
0.6 ms
0.9 ms
1.2 ms
1.5 ms
61
Figure 3.10: (a) The spatial distribution of shear stress σ
12
at time 1.5 ms and (b) snap-
shots of the slip profiles along the interface at various times for Case III with subshear
single-pulse rupture.
(a)
(b)
x
1
(m)
x
2
(m)
-2 -1 0 1 2
-2
-1
0
1
2
52
54
56
58
σ
12
-2 -1.5 -1 -0.5 0 0.5 1 1.5 2
0
0.1
0.2
0.3
0.4
0.5
0.6
Position along interface (m)
Slip (mm)
0.3 ms
0.6 ms
0.9 ms
1.2 ms
1.5 ms
62
Figure 3.11: (a) The spatial distribution of shear stress σ
12
at time 0.7 ms and (b) snap-
shots of the slip profiles along the interface at various times for Case IV with supershear
pulse-train rupture.
(a)
(b)
x
1
(m)
x
2
(m)
-2 -1 0 1 2
-2
-1
0
1
2
50
52
54
56
58
σ
12
-2 -1.5 -1 -0.5 0 0.5 1 1.5 2
0
0.2
0.4
0.6
0.8
1
1.2
1.4
1.6
1.8
Position along interface (m)
Slip (mm)
0.3 ms
0.5 ms
0.7 ms
0.9 ms
63
Figure 3.12: Simulations results with a reduced nominal dynamic coefficient of friction
µ
d
= 0.3 and the other parameters same as in Case II. (a) Profiles of frictional slip rate
(red) and shear traction (blue) along interface at time 0.5 ms. (b) Rupture front speeds V
R
vs. time. (c) Snapshots of the slip profile along the interface at various times. (d) Time
histories of various energy components with inset showing the approximate partition of
the total change of elastic energy between the kinetic energy W
kin
and heat energy W
heat
.
3.4 Discussion
In this chapter we examined numerically the energy partition and rupture properties
associated with different rupture modes, using two-dimensional in-plane finite element
calculations for a homogeneous isotropic elastic solid with an interface governed by a
RSD friction law proposed by Povirk and Needleman [1993]. Simulations based on vari-
ous initial shear stresses, rate dependence of the friction and nucleation parameters give
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7
0
0.5
1.0
1.5
2.0
2.5
3.0
Time (ms)
Rupture front speed V
R
(km/s)
β
β
α
c
R
2
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7
-3
-2
-1
0
1
2
3
Time (ms)
Energy (x10
6
J/m)
W
vol
el
Wapp
W
heat
W
int
el
W
kin
: 88%
: 12%
W
heat
W
kin
-2 -1 -0.5 0 0.5 1 1.5 2
30
40
50
54
60
Position along interface (m)
Shear stress (MPa)
-1.5
0
25
50
75
100
125
-2 -1.5 -1 -0.5 0 0.5 1 1.5 2
0
2
4
6
8
10
12
14
Position along the fault (m)
Slip (mm)
0.15 ms
0.3 ms
0.45 ms
0.6 ms
(a)
(b)
(c)
(d)
Slip rate (m/s)
64
rise to cases associated with supershear crack-like rupture, subshear crack-like rupture,
subshear single pulse and supershear train of pulses. The results add to the generality of
such rupture modes, and in particular to the mode of pulse train that appears to be transi-
tional between the crack and single-pulse modes of ruptures [Lapusta, 2005; Coker et al.,
2005; Shi and Ben-Zion, 2006].
A comparison between the results of the various simulated cases shows that increas-
ing the initial shear stress toward the failure stress can change the rupture speed from
subshear (Case II) to supershear (Case I), and increasing the rate dependence of the fric-
tion can change the rupture mode from a crack (Case II) to a pulse (Case III). These re-
sults are compatible with many previous studies of rupture on a planar frictional interface
in a homogeneous solid [e.g., Andrews, 1976; Cochard and Madariaga, 1996; Beeler and
Tullis, 1996; Zheng and Rice, 1998; Nielsen et al., 2000; Rosakis, 2002]. Our results also
indicate that the employed nucleation procedure can have significant effects on the rup-
ture mode. As pointed out by Ben-Zion [2006], relatively large nucleation zones with
relatively small stress drop favor crack-like ruptures, whereas relatively small nucleation
zones with relatively large stress drop favor pulse-like ruptures. In the simulations per-
formed in this work, the rupture mode switches from a single-pulse (Case III) to a train of
pulses (Case IV) when the width of the nucleation zone is doubled and the other condi-
tions are kept the same. Shi and Ben-Zion [2006] showed that the rupture speed along a
bimaterial frictional interface have a strong dependence on the strength of the nucleation
zone. Festa and Vilotte [2006] found in numerical studies with slip-weakening friction
65
that there is a trade-off between properties of the nucleation procedure and the initial
shear stress in determining the mode and speed of the dynamic ruptures.
The total release of elastic strain energy W
el
by different rupture modes over the same
propagation distance (Fig. 3.6) has the following order: supershear crack, subshear crack,
supershear train of pulses, subshear single pulse. The released elastic energy W
el
is parti-
tioned into radiated energy W
kin
and heat energy W
heat
. The ratio W
kin
/W
el
is referred to in
seismology as the seismic efficiency. In our study, with the same nominal static (0.6) and
dynamic (0.5) coefficients of friction, the seismic efficiency is ~6% for the supershear
crack-like rupture, ~3% for the subshear crack-like rupture, ~2% for the supershear train
of pulse, and only ~1% for the subshear single pulse. Reducing the nominal dynamic co-
efficient of friction to 0.3 and keeping the other parameters as in the subshear crack-like
rupture of Case II, increases the seismic efficiency to ~12% and produces a transition to
supershear rupture speed. Laboratory data associated with friction experiments in granite
[Lockner and Okubo, 1983] and estimates based on earthquake observations that cover a
wide range of event sizes and rock types [McGarr, 1999] indicate seismic efficiency val-
ues that are smaller than 6%. Kanamori et al. (1998) determined the source parameters of
the deep-focus 1994 great Bolivian earthquake (M
w
8.3) and suggested that the maximum
seismic efficiency was 3.6%. These estimates are compatible with the values obtained by
our simulated Cases I-IV.
The results indicate that crack-like ruptures are more efficient in terms of reducing
the stored elastic energy than pulse-like ruptures, and that supershear cracks are more ef-
ficient in this regard than subshear cracks. However, the subshear pulse occurs for a
66
lower background stress [e.g., Zhang and Rice, 1998] than the subshear, and especially
supershear cracks, and thus may be induced earlier than the subshear and supershear
cracks. In addition, the pulse mode of rupture is favored by heterogeneous stress condi-
tions [e.g., Das and Kostrov, 1988; Beroza and Mikumo, 1996; Day et al., 1998], while
crack-like ruptures require relatively homogeneous stress conditions [e.g., Ben-Zion and
Rice, 1997; Nielsen et al., 2000]. Earthquake faults are likely to be associated with highly
heterogeneous stress and they do not show evidence of the high localized dissipation that
may be expected for crack-like ruptures [e.g., Ben-Zion, 2001]. We also note that large
natural faults have bimaterial interfaces [e.g., McGuire and Ben-Zion, 2005] and that rup-
tures on bimaterial interfaces tend to propagate as subshear pulses [e.g., Weertman, 1980;
Ben-Zion, 2001; Ranjith and Rice, 2001; Shi and Ben-Zion, 2006; Ampuero and Ben-
Zion, 2006]. The above considerations suggest that the subshear pulse may be the most
common mode of rupture on earthquake faults, and that the supershear cracks may be the
least common. These expectations appear to be compatible with the available compila-
tions of earthquake data [e.g., Heaton, 1990; Ben-Zion, 2003; Mai, 2004; Manighetti et
al., 2005].
Continuing work on this topic may extend the current model to include bimaterial in-
terfaces [e.g., Andrews and Ben-Zion, 1997; Harris and Day, 1997; Shi and Ben-Zion
2006; Brietzke and Ben-Zion, 2006], plastic strain in the bulk [e.g., Andrews, 2005; Ben-
Zion and Shi, 2005] and spontaneous branching with creation of new surface area [e.g.,
Xu et al., 1997]. Simulations will incorporate ingredients that are characteristic of faults
at different evolutionary stage. The results can contribute to the improvement of our un-
67
derstanding of rupture processes in geologically relevant models and can be used to
check assumptions made in observational works.
68
Chapter 4
Slip Modes and Partitioning of Energy During
Dynamic Frictional Sliding Between Identical
Elastic-viscoplastic Solids
4.1 Introduction
Frictional sliding along an interface between two rapidly deforming solids is a basic
problem in mechanics encountered in a variety of contexts and at a variety of scales.
Rapid frictional sliding occurs on earthquake faults during propagation of earthquake
ruptures between two opposing fault blocks loaded by tectonic stresses. Complex fric-
tional sliding phenomena are also observed at much smaller scales [e.g., Persson, 2000].
Previous experimental and numerical studies, with a main focus on geophysical ap-
plications, have shown that frictional sliding on an interface between homogeneous elas-
tic solids can occur in various modes including a crack-like mode, a pulse-like mode and
a train-of-pulses mode, and also can propagate at various speeds including subshear and
supershear. In a crack-like mode, the slipping region continuously expands along the in-
terface and the slip duration at a given point in the sliding region is comparable to the
69
overall slip duration. On the other hand, in a pulse-like slip mode only a small portion of
the interface slides at a given time and the slip duration at a given point in the sliding re-
gion is much smaller than the overall sliding duration. Controlling factors on the selection
of sliding mode and propagation speed include: (i) the interface constitutive properties,
with a relatively strong dependence of friction on sliding speed favoring a pulse-like
mode and a relatively weak dependence favoring a crack-like mode [e.g., Beeler and Tul-
lis, 1996; Zheng and Rice, 1998; Nielsen et al., 2000; Ampuero and Ben-Zion, 2008]; (ii)
relatively high and smooth stress loading favors supershear propagation while relatively
low stress loading favors subshear rupture [e.g., Burridge, 19773; Andrews, 1976; Day,
1982; Rosakis, 2002; Xia et al., 2004; Shi et al., 2008]; and (iii) the strength, abruptness
and nature of the nucleation procedure used in earthquake modeling can have a strong
effect on both the sliding mode and the speed of propagation [e.g., Shi and Ben-Zion,
2006; Festa and Vilotte, 2006; Shi et al., 2008]. A contrast of elastic properties across the
fault can also affect the sliding mode and the propagation speed [e.g., Weertman, 1981;
Ben-Zion, 2001; Brietzke et al., 2009].
The analyses of frictional sliding mentioned above assume that the bulk material is
elastic. However, the elastically predicted stress concentration near the sliding interface
may induce the bulk material to respond inelastically [Poliakov et al., 2002; Rice et al.,
2005]. The generation of inelastic strain may affect the sliding mode and the speed of
propagation as well as the partitioning of energy. Several numerical studies have been
carried out in the geophysics literature on the effects of off-fault inelastic response during
crack-like dynamic ruptures of earthquake faults. These studies focused on the damage
70
pattern in the bulk and treated off-fault inelastic response in various forms including for-
mation of tensile cracks [e.g., Yamashita, 2000; Dalguer et al., 2003], shear branches
[e.g., Ando and Yamashita, 2007] or generation of plastic strain using a Coulomb yield-
ing criterion [e.g., Andrews, 2005; Templeton and Rice, 2008]. Off-fault Coulomb yield-
ing was also used in dynamic simulations of a pulse-like rupture on a frictional interface
between two solids with different properties [Ben-Zion and Shi, 2005].
Here, we carry out dynamic analyses of mode-II frictional sliding along an interface
between two identical solids under plane strain conditions initiated by impact loading.
The calculations employ the same configuration analyzed in Coker et al. [2005] but here
plane strain conditions are considered. Two identical rectangular regions are held to-
gether by an initial compressive stress. Shear loading is imposed by edge impact near the
interface. Several loading conditions are considered that give rise for purely elastic be-
havior to frictional slip in a crack-like mode, in a pulse-like mode and in a train-of-pulses
mode. The materials are modeled as isotropically hardening and rate dependent. Friction
is characterized by a rate- and state-dependent relation [Dieterich, 1979; Ruina, 1983;
Rice and Ruina, 1983; Linker and Dieterich, 1992; Povirk and Needleman, 1993], includ-
ing the Prakash-Clifton [Prakash and Clifton, 1993; Prakash, 1998] normal stress de-
pendence. Our investigation focuses on the effects of plasticity on the mode and speed of
frictional sliding along the interface and on the partitioning of energy. We also carry out
parameter studies exploring the effects of variations in interface elastic stiffness and im-
pact velocity.
71
4.2 Problem Formulation
4.2.1 Boundary Value Problem Formulation
The finite element formulation is the same as in Chapter 3. In Cartesian tensor nota-
tion, the principle of virtual work is written as
€
σ
ij
:δε
ij
dV
V
∫
− T
i
δΔu
i
dS
S
int
∫
= T
i
δΔu
i
dS
S
ext
∫
− ρ
∂
2
u
i
∂t
2
δu
i
dV
V
∫
(4.1)
where t is time, σ
ij
are the stress tensor components, u
i
are the displacement vector com-
ponents, ε
ij
are the components of the infinitesimal strain tensor, ∆ denotes a jump across
an interface, and V, S
ext
and S
int
are the volume, external surface area and internal surface
area, respectively. The density of the material is ρ and T
i
= σ
ij
n
j
with n
j
being the surface
normal components.
Figure 4.1: Geometry of the configuration analyzed and time function of the impact load-
ing. The dashed line denote interface S
int
(x
2
= 0) where frictional sliding occur.
Σ
0
x
2
2l
w
b
V(t)
t
V(t)
V
imp
t
r
S
int
t
p
t
s
Σ
0
x
1
72
Computations are carried out for the specimen geometry shown in Fig. 4.1 with l = 75
mm and w = 132 mm. Plane strain conditions are assumed and a Cartesian coordinate
system is used with the x
1
−x
2
plane as the plane of deformation. For t ≤ 0, the body is
subject to a uniform initial compressive stress Σ
0
and this stressed state is used as a refer-
ence for deformation measurements. A normal velocity V(t) is prescribed on the edge x
1
=
0, along the interval −b ≤ x
2
≤ 0, with b = 25 mm, where the shear traction vanishes, so
that
€
u
1
=− V ξ
( )
dξ
0
t
∫
and T
2
= 0 on x
1
= 0 and –b ≤ x
2
≤ 0. (4.2)
In Eq. (4.2)
€
V t ( ) =
V
imp
t /t
r
, for 0≤ t < t
r
V
imp
, for t
r
≤ t≤ t
p
+ t
r
( )
V
imp
1− t−t
p
−t
r
( ) [ ]
t
s
, for t
p
+ t
r
( )
< t < t
s
+ t
p
+ t
r
( )
0, for t≥ t
s
+ t
p
+ t
r
( )
(4.3)
where, t
r
is the rise time, t
p
is the pulse width and t
s
is the step-down time. On the remain-
ing external surfaces, the tractions vanish, i.e. T
i
= 0.
The finite element discretization and integration scheme are the same as in Chapter 3.
The mesh employed is the same as in Coker et al. [2005] which consists of a fine uniform
region near the interface close to the side of impact connected to a coarse outer region by
a transitional region with gradually varying mesh spacing (Fig. 4.2). The domain ana-
lyzed has 18,320 quadrilateral elements and 74,360 degrees of freedom. The uniform
fine-mesh region is 60 mm × 4 mm and is discretized into 4000 quadrilaterals, each of
which is 0.3 mm × 0.2 mm. Along the interface S
int
, four-node elements are used to give a
73
linear dependence of the displacement jump on x
1
. The bulk finite elements are constant
strain elements and therefore, effectively, one-point integration is used to form the ele-
ment stiffness matrices. The interface element stiffness matrices are integrated using
four-point Gaussian quadrature.
Figure 4.2: Close-up view of the mesh near the transition x
1
= 60 mm between the region
of fine uniform mesh and the region of graduated mesh along S
int
.
4.2.2 Bulk elastic-viscoplastic constitutive relation
The bulk material is characterized as an elastic-viscoplastic isotroically hardening
solid. The total strain rate,
€
˙ ε
ij
e
, and a plastic part,
€
˙ ε
ij
p
, such that
€
˙ ε
ij
= ˙ ε
ij
e
+ ˙ ε
ij
p
(4.4)
with the overdot denoting partial differentiation with respect to time and
€
˙ ε
ij
e
=
1+ν
E
˙ σ
ij
−
ν
E
δ
ij
˙ σ
kk
,
€
˙ ε
ij
p
=
3
2
˙
ε
σ
e
′ σ
ij
. (4.5)
x
1
40 50 60 80
−20
−10
0
10
(mm)
x
2
(mm)
70
20
74
In Eq. (4.5), δ
ij
is the Kronecker delta, E is the Young’s modulus, ν is the Poisson’s
ratio,
€
˙
ε is the effective plastic strain rate and σ
e
is the Mises effective stress defined as
€
˙
ε =
2
3
˙ ε
ij
p
˙ ε
ij
p
(4.6)
€
σ
e
=
3
2
′ σ
ij
′ σ
ij
(4.7)
where
€
′ σ
ij
is the deviatoric stress given by
€
′ σ
ij
=σ
ij
−
1
3
σ
kk
δ
ij
(4.8)
and
€
˙
ε = ˙ ε
0
σ gε
( ) [ ]
1m
,
€
gε
( )
=σ
0
1+ε ε
0
( )
N
.
(4.9)
In Eq. (4.9), σ
0
and ε
0
are the reference stress and strain related by σ
0
=
Eε
0
and
€
ε =
˙
ε dt
∫
is the effective plastic strain. The function
€
gε
( )
specifies the strain hardening
characteristics with σ
e
=
€
gε
( )
representing the effective stress versus effective strain re-
sponse in a tensile test carried out at
€
˙
ε equal to the reference strain rate
€
˙ ε
0
, N is the strain
hardening exponent and m is the strain-rate hardening exponent with m 0 correspond-
ing to the rate-independent limit. A rate tangent modulus method is used for the visco-
plastic stress update [Peirce et al., 1984].
Since the material response is viscoplastic, the speeds of dilational, shear and
Rayleigh waves are still given by the elastic expressions [e.g., Freund, 1998]
€
c
l
=
E 1−ν ( )
ρ 1+ν ( ) 1−2ν ( )
,
€
c
s
=
E
2ρ 1+ν ( )
,
€
c
R
=c
s
0.862 +1.14ν
1+ν
. (4.10)
75
As in Coker et al. [2005], the elasticity parameters are chosen to represent Homalite-
100: E = 5.3 GPa, ν = 0.35 and ρ = 1246 kg/m
3
. The corresponding wave speeds given by
Eq. (4.10) are c
l
= 2613 m/s, c
s
= 1255 m/s and c
R
= 1174 m/s.
4.2.3 Interface constitutive relations
As in Chapter 3, the interface S
int
along x
2
= 0 is characterized by an elastic and a
rate-dependent frictional relation [Dieterich, 1979; Ruina, 1983; Rice and Ruina, 1983;
Linker and Dieterich, 1992; Povirk and Needleman, 1993] specified by
€
˙
T
n
=−C
n
Δ˙ u
n
, (4.11)
€
˙
T
s
=−C
s
Δ˙ u
s
− sgn T
s
( )Δ˙ u
slip [ ]
. (4.12)
Here, the subscripts s and n denote quantities tangential and normal to the interface,
respectively, and C
s
and C
n
are the interface stiffness parameters. Eqs (4.11) and (4.12)
apply if the surfaces are in contact; if they separate, then T
n
= T
s
= 0.
The jump of the tangential particle velocity across the interface,
€
Δ˙ u
s
, consists of an
elastic part
€
Δ˙ u
s
el
and a frictional part
€
Δ˙ u
slip
such that
€
Δ˙ u
s
=Δ˙ u
s
el
+Δ˙ u
slip
(4.13)
where
€
Δ˙ u
s
el
=−
˙
T
s
C
s
.
We use the phenomenological expression for
€
Δ˙ u
slip
proposed by Povirk and Needle-
man [1993] with a Prakash-Clifton regularization for normal stress dependence. In Eq.
(4.13),
€
Δ˙ u
slip
= 0 if T
n
≤ 0, while for T
n
> 0, which corresponds to the interface being un-
der compression,
€
Δ˙ u
slip
is given by
76
€
Δ˙ u
slip
= V
0
T
s
θ
PC
g θ
0
( )
r
−1
for
€
T
s
θ
PC
gθ
0
( )
>1, (4.14)
and
€
Δ˙ u
slip
= 0 otherwise. In Eq. (4.14)
€
gθ
0
( )
=
µ
d
+ µ
s
− µ
d
( )
exp − L
0
θ
0
( )
V
1
( )
p
[ ]
1+ L
0
θ
0
( )
V
0
( )
1 r
(4.15)
where the evolution of the state variable θ
0
is given by
€
˙
θ
0
= B 1−
θ
0
Δ˙ u
slip
L
0
(4.16)
and the evolution of the Prakash-Clifton internal variable θ
PC
is given by
€
˙
θ
PC
=−
1
L
PC
˙
θ
PC
− T
n ( )
Δ˙ u
slip
. (4.17)
In Eqs (4.12)-(4.17) µ
s
and µ
d
are the nominal static and dynamic coefficients of fric-
tion, L
0
is the characteristic slip distance over which the friction coefficient evolves to a
steady-state value, L
PC
is the characteristic slip distance over which the normal stress
evolves to a steady-state value, V
0
and V
1
are reference slip velocities. All the calculations
in this chapter employ the parameters values µ
s
= 0.6, µ
s
= 0.5, V
0
= 100 m/s, V
1
= 26
m/s, L
0
= L
PC
= 0.2 mm, r = 5, p = 1.2 and B = 4.6. Unless otherwise specified, C
n
= 300
GPa/m, C
s
= 100 GPa/m. As noted in the previous Chapter, only a subset of those values
was found in previous studies to characterize experimental observations of Homalite slid-
ing on Homalite [e.g., Rosakis et al, 2000; Rosakis, 2002; Samudrala et al., 2002]. A re-
stricted form of the Prakash-Clifton response with only one state variable is employed.
77
The interface shear traction T
s
is updated using a rate tangent method the formula for
which is given in the Appendix (Section 4.5).
4.2.4 Energy Partitioning
One focus of this study is the effect of plasticity on energy partitioning which is of
particular interest in earthquake dynamics [e.g., Husseini, 1977; Shi et al. 2008]. Follow-
ing Chapter 3, the energy balance equation can be written as
€
W
boun
−W
vol
el
−W
vol
pl
−W
kin
+φ
el
+φ
slip
=0 (4.18)
where the various energy components are defined as
(i) Energy input at the external boundaries
€
W
boun
= T
i
˙ u
i
dS
S
ext
∫
[ ]
dt
0
t
∫
(4.19)
(ii) Change of elastic strain energy and plastic dissipation in the bulk
€
W
vol
el
= σ
ij
˙ ε
ij
e
dV
V
∫ [ ]
dt
0
t
∫
,
€
W
vol
pl
= σ
ij
˙ ε
ij
p
dV
V
∫ [ ]
dt
0
t
∫
(4.20)
(iii) Kinetic energy (or radiated energy)
€
W
kin
=
1
2
ρ˙ u
i
˙ u
i
dV
V
∫
(4.21)
(iv) Change of interface energy consisting a frictional part and an elastic part
€
φ
slip
= T
i
Δ˙ u
slip
dS
S
int
∫
[ ]
dt
0
t
∫
(4.22)
€
φ
el
= T
s
Δ˙ u
s
− sgn T
s
( )Δ˙ u
slip [ ]
dS
S
int
∫
+ T
n
Δ˙ u
n
dS
S
int
∫
[ ]
dt
0
t
∫
(4.23)
78
4.3 Numerical Results
Calculations were carried out for three combinations of initial compressive stress Σ
0
and impact velocity V
imp
considered in Coker et al. [2005]
Case I: A crack-like mode – Σ
0
= 6 MPa, Σ
e
= 5.27 MPa, V
imp
= 2 m/s.
Case II: A pulse-like mode – Σ
0
= 30 MPa, Σ
e
= 26.4 MPa, V
imp
= 2 m/s.
Case III: A train-of-pulses mode – Σ
0
= 10 MPa, Σ
e
= 8.79 MPa, V
imp
= 20 m/s.
where Σ
e
is calculated from Σ
0
using Eq. (4.7).
Our aim is to explore the effect of bulk material inelasticity on the mode of frictional
sliding and on the partitioning of energy rather than to describe the response of a particu-
lar material. For each case we present results from calculations with elastic material be-
havior (hereafter referred to as “elastic calculations”) and results form calculations where
the bulk material reference strength σ
0
is taken to be slightly greater than Σ
e
, the initial
effective stress associated with the initial stress state (hereafter referred to as “plastic cal-
culations”). We employ a simple isotropic hardening material characterization and, since
the initial stress state differs in each case, the material description also varies from case to
case. The term used to characterize each case is the frictional slip mode obtained in the
elastic calculations.
The strain hardening exponent N is fixed at 0.1 and the reference strain rate is
€
˙ ε
0
= 1
sec
−1
. In Cases I and III, the strain-rate hardening exponent m is taken to be 0.01 giving a
nearly rate-independent response, while in Case II, for numerical stability reasons, m is
set to 0.05 giving a stronger rate dependence. The calculations are carried out using an
initial time step of 0.05 × min(ds)/c
l
in Cases I and III and 0.002 × min(ds)/c
l
in Case II,
79
where min(ds) denotes the shortest edge length of the triangular elements. An adaptive
time step scheme is employed that limits the maximum value of
€
Δ˙ u
slip
dt in a single time
step to 0.01L
0
.
In all cases an elastic wave propagates into the bulk material and along interface S
int
with the dilational wave speed c
l
at the initiation of impact. The loading wave induces an
elastic displacement jump,
€
Δ˙ u
s
el
, along interface S
int
and a corresponding interfacial shear
traction T
s
. After some time, which depends on the loading condition and on the proper-
ties of the bulk material and the interface, frictional sliding initiates and propagates along
the pre-loaded interface S
int
. The calculations are terminated before the frictional slip
front reaches the end of the uniform fine-mesh region, i.e., x
1
= 60 mm. We calculate en-
ergy changes both in the bulk and on the interface from the state at t = 0.
4.3.1 Case I: a crack-like mode
(a) (b)
Figure 4.3: Along-interface profiles of
€
Δ˙ u
slip
and T
s
/T
n
in Case I at 20, 25 and 30 µs for
(a) the elastic calculation and (b) the plastic calculation.
With Σ
0
= 6 MPa and V
imp
= 2 m/s, a crack-like slip mode is obtained in both the elas-
tic calculation and the plastic calculation with σ
0
= 5.3 MPa (slightly above Σ
e
). Fig. 4.3
0 10 20 30 40 50 60 70
0
4
8
x
1
(mm)
Δu
slip
(m/s)
20 μs
25 μs
30 μs
6
2
0
0.2
0.4
0.6
T
s
/
T
n
0 10 20 30 40 50 60 70
0
2
4
6
8
10
12
Δu
slip
(m/s)
x
1
(mm)
20 μs
25 μs
30 μs
0
0.2
0.4
0.6
T
s
/
T
n
80
plots the frictional slip rate,
€
Δ˙ u
slip
, along the interface together with the apparent coeffi-
cient of friction, µ
app
= T
s
/T
n
, at different times, for both elastic and plastic calculations.
The slip rate in the plastic calculation is greater than in the elastic calculation and there is
a pronounced peak in slip rate near the frictional slip front. In addition, in the plastic cal-
culation, there is a small region along S
int
near the edge of impact where
€
Δ˙ u
slip
= 0 due to
a loss of contact so that T
n
= T
s
= 0 and µ
app
is undefined.
In the elastic calculation, frictional sliding along S
int
nucleates from the boundary near
the impact edge (x
1
= 0) at t = 15.72 µs when T
s
increases to 3.6 MPa and T
n
decreases to
5.7 MPa so that µ
app
= 0.63 there. In the plastic calculation, frictional sliding nucleates
earlier, at t = 12.39 µs, with T
n
≈ T
s
= 3.6 MPa (µ
app
≈ 1) also at the impact end of the
interface. In both calculations,
€
Δ˙ u
slip
increases sharply after nucleation while
€
Δ˙ u
s
el
de-
creases sharply to slightly below zero and remains small over most of the region of fric-
tional sliding. Both T
s
and T
n
in the plastic calculation are smaller everywhere along S
int
than in the elastic calculation and vanish near the impact end.
Fig. 4.4 shows curves of frictional slip propagation speeds V
slip
for both calculations.
Similar to the procedure in Chapter 3, the frictional slip propagation speed is calculated
by first determining the position of a frictional slip front, defined as the furthest distance
from the impact end where
€
Δ˙ u
slip
exceeds a specified threshold value. A seven-point pro-
gressive least-squares fit to the slip front position as a function of time is then used to
calculate V
slip
. In all the calculations in this Chapter, the threshold value of
€
Δ˙ u
slip
is cho-
sen to be 5 × 10
−4
m/s. Calculations using other threshold values of
€
Δ˙ u
slip
show little dif-
81
ference in the calculated frictional slip propagation speeds. This insensitivity to the
choice of threshold value of
€
Δ˙ u
slip
is due to the sharpness of the frictional slip front as
shown in Fig. 4.3. For the elastic calculation, the calculated frictional slip propagation
speed initially significantly exceeds the elastic longitudinal wave speed but then de-
creases quickly as slip propagates and attains a quasi-steady value close to c
l
. On the
other hand, the calculated frictional slip propagation speed in the plastic calculation is
initially less than c
l
but increases quickly to a value close to but slightly less than that of
the elastic calculation. As noted in Coker et al. [2005], the initiation process of frictional
sliding under impact loading can involve coalescence of small sliding regions. Therefore
the calculated apparent propagation speed of frictional slip front during the early stages
can exceed c
l
.
Figure 4.4: Propagation speed of the frictional sliding front, V
slip
, versus time in Case I
for the elastic and plastic calculations.
Fig. 4.5 shows time histories of various energies including the change of elastic strain
energy in the bulk
€
W
vol
el
, the plastic dissipation in the bulk
€
W
vol
pl
, the kinetic energy
€
W
kin
,
5 10 15 20 25 30 35 40
1
2
3
4
5
6
7
Time (μs)
elastic
plastic
c
l
V
slip
(km/s)
82
and the frictional dissipation along the interface
€
φ
slip
, for both elastic and plastic calcula-
tions. In the elastic calculation (Fig. 4.5a), during the early stages of loading before the
initiation of frictional slip, the energy input is divided between
€
W
vol
el
and
€
W
kin
. Once fric-
tional slip initiates, the energy dissipation from frictional sliding increases and eventually
exceeds both
€
W
vol
el
and
€
W
kin
. The energy balance involving all the energy terms in Eq.
(4.18) was checked in this and all other calculations and the mismatch is less than 0.5%
of the bulk energy (
€
W
vol
el
+
€
W
vol
pl
). Fig. 4.5(b) shows that in the plastic calculation the bulk
elastic energy initially decreases due to plastic relaxation, and then increases somewhat,
but the bulk energy is mainly plastic dissipation.
(a) (b)
Figure 4.5: Time histories of various energy components in Case I for (a) the elastic cal-
culation and (b) the plastic calculation.
Contours of Mises effective stress σ
e
, defined in Eq. (4.8), are shown in Fig. 4.6(a)
for the elastic calculation and in Fig. 4.6(b) for the plastic calculation. Plasticity limits the
effective stress levels that can be attained so that the scale in Fig. 4.6(b) is lower than that
in Fig. 4.6(a). In both calculations, the highest values of σ
e
occur for x
2
> 0 which is the
half of the body above the impact surface. Stress relaxation occurs at the interface near
0 5 10 15 20 25 30 35 40
−2
0
2
4
6
8
10
12
14
Time (μs)
Energy (J/m)
W
vol
φ
slip
W
vol
W
kin
el
pl
0 5 10 15 20 25 30 35 40
−2
0
2
4
6
8
10
12
14
16
18
Time (μs)
Energy (J/m)
W
vol
φ
slip
W
vol
W
kin
el
pl
83
the impact surface in Fig. 4.6(b), which leads to the arrest of frictional sliding in a small
region. Although contours of plastic strain rate and plastic strain for Case I are not
shown, the maximum plastic strain is about 0.01 and the maximum plastic strain rate is
about 2 × 10
3
sec
−1
.
(a)
(b)
Figure 4.6: Contours of Mises effective stress in Case I for (a) the elastic calculation at t
= 34.6 µs and (b) the plastic calculation at t = 34.6 µs. The frictional sliding front is at x
1
= 58.0 mm in (a) and at x
1
= 55.5 mm in (b).
x
1
0 20 40 60 80
−40
−20
0
20
(mm)
x
2
(mm)
14
12
10
8
6
4
2
σ
e
(MPa)
6
5
4
3
2
1
x
1
0 20 40 60 80
−40
−20
0
20
(mm)
x
2
(mm)
σ
e
(MPa)
84
4.3.2 Case II: a pulse-like mode
In Case II, the initial stress is increased to Σ
0
= 30 MPa with the impact velocity re-
maining at V
imp
= 2 m/s. For numerical stability reasons, the rate sensitivity parameter in
the plastic calculation with σ
0
= 27 MPa is taken to be m = 0.05, thus increasing the mate-
rial strain rate sensitivity.
(a) (b)
Figure 4.7: Along-interface profiles of
€
Δ˙ u
slip
and T
s
/T
n
in Case II at 30, 35 and 40 µs for
(a) the elastic calculation and (b) the plastic calculation.
In the elastic calculation frictional sliding occurs in a pulse-like mode as shown in
Fig. 4.7(a). Instead of nucleating at the boundary as in Case I, frictional sliding in the
elastic calculation initiates as a small slip patch at t = 29.55 µs at about x
1
= 5.7 mm when
T
s
reaches a peak value of 18 MPa there. The slip patch expands in both directions with
an asymmetric slip-rate profile in a shape that has a gentle slope away from the loading
edge and a steep slope towards the loading edge. The slip-rate amplitude near the steep-
slope end of the frictional slip patch increases giving the slip patch a pulse-like shape as it
approaches the impact edge while the gentle-slope end of the slip patch maintains its
shape. At t = 30.73 µs, the pulse-like slip patch meets the boundary with a peak frictional
slip rate of 1325 m/s and its single peak starts to split into two peaks creating two pulses
0 10 20 30 40 50 60 70
0
200
400
600
800
1000
Δu
slip
(m/s)
0
0.2
0.4
0.6
T
s
/
T
n
x
1
(mm)
30 μs
35 μs
40 μs
slip patch
0 10 20 30 40 50 60 70
0
200
400
600
800
1000
1200
1400
1600
Δu
slip
(m/s)
0
0.2
0.4
0.6
T
s
/
T
n
30 μs
35 μs
40 μs
x
1
(mm)
85
going in opposite directions: one exits the impact end of the boundary and the other
propagates away from that boundary with a peak frictional slip rate of about 600 m/s. The
shear traction (and the shear traction rate) has a local minimum at the location of the slip
pulse. As the frictional slip pulse propagates, the local minimum of shear traction moves
with the pulse and the amplitude of a second peak of shear traction that follows behind
the pulse gradually increases to nucleate another frictional slip patch that leads to a sec-
ond pulse. This process repeats to generate more pulses. The nucleation location of the
slip patch moves away from the impact end: the second, third and fourth slip patches nu-
cleate at about x
1
= 6.5, 8.5 and 12.5 mm. A developing slip patch can be seen in the slip-
rate profile at 30 µs in Fig. 4.7(a).
As shown in Fig. 4.7(b), the pulse-like mode of propagation is preserved in the plastic
calculation and the initiation of frictional occurs in a similar fashion to that in the elastic
calculation. A small frictional slip patch nucleates around x
1
= 2.5 mm at t = 28.32 µs
when T
s
reaches a peak value of 18 MPa there. Since the first slip patch is closer to the
impact edge, it reaches the boundary only 0.63 µs after nucleation at t = 28.95 µs, when
the peak slip rate exceeds 10000 m/s. The second slip patch nucleates at about the same
location as the first slip patch when the second peak of T
s
reaches 16 MPa and µ
app
reaches 0.65 there. The third and ensuing slip patches, however, nucleate at the boundary
due to the rapid decrease of T
n
there which gives rise to a large value of µ
app
. The time
interval between the nucleation of slip patches decreases leading to a smaller spacing be-
tween pulses. Since the initial and subsequent slip patches have less time to grow than in
the elastic calculation (because they nucleate closer to the boundary), the slip pulse has a
86
smaller width than in the elastic calculation. By the time the first frictional slip pulse
reaches the end of the fine-mesh region, there are eight fully developed pulses and one
developing pulse in the plastic calculation, while in the corresponding elastic calculation
there are four fully developed pulses and one developing pulse.
For both elastic and plastic calculations the shear traction rate
€
˙
T
s
is negative at the lo-
cation of the frictional slip pulses and therefore so is
€
Δ˙ u
s
el
. In fact, as also seen by Coker
et al. [2005], the magnitude of
€
Δ˙ u
s
el
is such that the velocity jump across the interface,
€
Δ˙ u
s
, exhibits a rather smooth variation along the interface despite the pulse nature of
€
Δ˙ u
slip
.
The numerical issue leading to the use of m = 0.05 is due to high pulse frictional slip
rates and corresponding high bulk plastic strain rates which necessitate extremely small
time steps for numerical stability (with the algorithms used here). The plastic strain rate
increases with decreasing m. For example, with m = 0.05 as in the calculation here the
peak plastic strain rate is about 10
4
sec
−1
, while with m = 0.01 the peak plastic strain rate
can reach 10
8
sec
−1
.
For both elastic and plastic calculations the calculated frictional slip propagation
speed significantly exceeds c
l
over the entire time interval calculated (Fig. 4.8).
Fig. 4.9 shows the rapid rise in the energy dissipated in slip,
€
φ
slip
, once frictional slip
propagation initiates. As for the elastic calculation in Case I, the values of
€
W
vol
el
and
€
W
vol
pl
remain close over the entire range of the computation. The absolute values of the various
energy contributions are not directly comparable to those in Case I because the initial
stress states differ.
87
Figure 4.8: Propagation speed of the frictional sliding front, V
slip
, versus time in Case II
for the elastic and plastic calculations.
(a) (b)
Figure 4.9: Time histories of various energy components in Case II for (a) the elastic cal-
culation and (b) the plastic calculation.
Comparing the energy partitioning in Fig. 4.9 shows that in the plastic calculation the
bulk stress working is mainly plastic dissipation rather than stored elastic energy. The
time evolution of kinetic energy is not much changed. Bulk material plasticity results in
greater frictional dissipation along the interface at a given time but this is mainly a shift
25 30 35 40 45
1
2
3
4
5
6
7
8
9
10
Time (μs)
elastic
plastic
V
slip
(km/s)
c
l
0 5 10 15 20 25 30 35 40 45
−20
0
20
40
60
80
100
120
140
160
180
Time (μs)
Energy (J/m)
W
vol
φ
slip
W
vol
W
kin
el
pl
0 5 10 15 20 25 30 35 40 45
−20
0
20
40
60
80
100
120
140
160
180
Time (μs)
Energy (J/m)
W
vol
φ
slip
W
vol
W
kin
el
pl
88
due to the earlier nucleation time of frictional slip in the plastic calculation than in the
elastic calculation.
(a)
(b)
Figure 4.10: Contours of Mises effective stress in Case II for (a) the elastic calculation at
t = 37.6 µs and (b) the plastic calculation at t = 37.8 µs. The frictional sliding front is at x
1
= 42.2 mm in (a) and at x
1
= 42.4 mm in (b).
The distributions of Mises effective stress σ
e
in Fig. 4.10 reflect the pulse nature of
frictional sliding for both elastic and plastic calculations. As in Case I, the effective stress
48
46
42
40
35
30
25
x
1
0 20 40 60 80
−20
0
20
(mm)
x
2
(mm)
−40
σ
e
(MPa)
38
36
34
32
30
28
26
x
1
0 20 40 60 80
−40
−20
0
20
(mm)
x
2
(mm)
σ
e
(MPa)
89
levels are somewhat reduced in the plastic calculation and the high values of σ
e
mainly
exist in the region x
2
> 0.
4.3.3 Case III: a train of pulses
In Case III, the initial compressive stress Σ
0
is 10 MPa but the impact velocity is 20
m/s, an order of magnitude larger than in Cases I and II. The bulk material reference
strength for the plastic calculation σ
0
set to be 10 MPa.
(a) (b)
Figure 4.11: Along-interface profiles of
€
Δ˙ u
slip
and T
s
/T
n
in Case III (a) at 10, 20 and 29.4
µs for the elastic calculation and (b) at 10, 20, 30 µs for the plastic calculation.
In the elastic calculation, frictional sliding nucleates near the boundary at t = 7.06 µs
and starts to propagate in a crack-like fashion. The value of µ
app
at the boundary continue
to increase until t = 11.42 µs when it reaches a peak value of 0.54. Then T
s
abruptly de-
creases while T
n
remains relatively constant leading to a sharp decrease of µ
app
to about
0.5, which results in
€
Δ˙ u
slip
vanishing and leaving the initial crack-like region detached
from the following slip. At the boundary, T
n
monotonically decreases while T
s
oscillates
with an overall decreasing amplitude, which leads to µ
app
oscillating between 0.46 and
0.62 generating a train of pulses as shown in Fig. 4.11(a). Unlike the situation in Case I
0 10 20 30 40 50 60 70
0
100
200
300
400
0
0.2
0.4
0.6
T
s
/
T
n
Δu
slip
(m/s)
x
1
(mm)
10 μs
20 μs
29.4 μs
0 10 20 30 40 50 60 70
50
100
150
200
250
0
Δu
slip
(m/s)
10 μs
20 μs
30 μs
0
0.2
0.4
0.6
T
s
/
T
n
x
1
(mm)
90
and Case II, T
n
remains above its initial value along the entire interface S
int
throughout the
calculation.
In the plastic calculation, frictional sliding initiates from the boundary earlier at t =
6.3 µs. Also the oscillation of T
s
occurs at a much shorter interval after initiation with an
amplitude not large enough to bring
€
Δ˙ u
slip
to zero. Therefore the frictional slip mode re-
mains crack-like overall but with strong oscillations near the leading front. A region with
vanishing T
s
and T
n
occurs near the impact end of the interface at t = 13.6 µs. As slip con-
tinues to propagate, the oscillations of T
s
increase eventually transforming the oscillatory
portion of the crack into a train of pulses as shown in Fig. 4.11(b).
The peak frictional slip rate in a pulse in both calculations is about half that for the
pulse-like propagation in Case II (Section 4.3.2). In the plastic calculation, the peak fric-
tional slip rate at t = 30 µs in the crack-like slip region is a factor of 4-5 times smaller
than in the region with a train of pulses, and the pulses near the leading frictional slip
front are more closely spaced than in the corresponding elastic calculation.
Compared to Fig. 4.4 in Case I, the initial evolution of the frictional slip propagation
speeds in Case III for both elastic and plastic calculations shown in Fig. 4.12 are similar.
In Case III the propagation speeds are more oscillatory and slightly slower with mean
values of about 2% to 4% below c
l
.
For the elastic calculation shown in Fig. 4.13(a), as in the previous elastic calcula-
tions, the bulk elastic energy
€
W
vol
el
and the kinetic energy
€
W
kin
are almost equal. How-
ever, for the train-of-pulses mode, the frictional dissipation
€
φ
slip
remains much less than
either
€
W
vol
el
or
€
W
kin
over the entire time range of the computation.
91
Figure 4.12: Propagation speed of the frictional sliding front, V
slip
, versus time in Case III
for the elastic and plastic calculations.
(a) (b)
Figure 4.13: Time histories of various energy components in Case III for (a) the elastic
calculation and (b) the plastic calculation.
For the plastic calculation shown in Fig. 4.13(b), the bulk elastic and kinetic energies
dominate the energy distribution during the early stages of propagation (t ≤ 10 µs) but in
the latter stages, e.g., t > 25 µs, when the crack-like region has become more prominent,
the bulk plastic dissipation and frictional dissipation on the interface dominate. At t = 30
0 5 10 15 20 25 30
1.5
2
2.5
3
3.5
4
Time (μs)
V
slip
(km/s)
c
l
elastic
plastic
0 5 10 15 20 25 30
0
50
100
150
200
250
300
350
400
450
Time (μs)
Energy (J/m)
W
vol
φ
slip
W
vol
W
kin
el
pl
0 5 10 15 20 25 30 35
0
50
100
150
200
250
300
350
400
450
Time (μs)
Energy (J/m)
W
vol
φ
slip
W
vol
W
kin
el
pl
92
µs, the frictional dissipation is about a factor of two greater than it is for the correspond-
ing elastic calculation and the kinetic energy is similarly reduced.
(a)
(b)
Figure 4.14: Contours of Mises effective stress in Case III for (a) the elastic calculation at
t = 29.4 µs and (b) the plastic calculation at t = 30 µs. The frictional sliding front is at x
1
= 59.1 mm in (a) and at x
1
= 58.2 mm in (b).
The distribution of Mises effective stress in Fig. 4.14(a) for the elastic calculation re-
flects the train-of-pulses propagation mode of frictional slip. In contrast, for the plastic
42
41
40
35
30
25
20
15
10
5
x
1
0 20 40 60 80
−40
−20
0
20
(mm)
x
2
(mm)
σ
e
(MPa)
14
12
10
8
6
4
2
x
1
0 20 40 60 80
−40
−20
0
20
(mm)
x
2
(mm)
σ
e
(MPa)
93
calculation, the oscillations that are indicative of pulses are smoothed out and the distri-
bution resembles that of the crack-like mode as in Fig. 4.6(b). Larger plastic strains de-
velop in Case III than in Case II (values of effective plastic strain
€
ε greater than 0.10)
and plastic straining also occurs over a larger region acting to smooth the oscillations.
4.3.4 Parameter studies
4.3.4.1 Effect of interface elasticity
(a) (b)
Figure 4.15: Along-interface profiles of
€
Δ˙ u
slip
with the frictional slip front position at 50
mm in (a) the elastic calculations and (b) the plastic calculations. V
imp
= 20 m/s and C
s
=
100, 200 and 400 GPa/m. The time for each profile is indicated.
In order to explore the effect of interface elasticity, we repeat the calculations in Case
III using various values for the interface stiffness parameters C
s
and C
n
. Fig. 4.15 shows
the slip-rate profiles from elastic calculations and from plastic calculations with σ
0
= 10
MPa using C
s
= 100 (the value used in the calculations in Section 4.3.3), 200 and 400
GPa/m with all other parameters the same as in Section 4.3.3. For the elastic calculations
shown in Fig. 4.15(a), frictional slip propagates as a train of pulses with C
s
= 100 GPa/m.
When C
s
is increased to 200 GPa/m, the frictional slip mode changes to a crack-like
0 10 20 30 40 50 60
0
50
100
150
200
250
300
350
Cs = 100 GPa, t = 25.72 μs
Cs = 200 GPa, t = 23.34 μs
Cs = 400 GPa, t = 21.74 μs
Δu
slip
(m/s)
x
1
(mm)
0 10 20 30 40 50 60
0
50
100
150
200
250
Cs = 100 GPa, t = 26.57 μs
Δu
slip
(m/s)
x
1
(mm)
Cs = 200 GPa, t = 21.95 μs
Cs = 400 GPa, t = 23.70 μs
94
mode with notable oscillations in the slip-rate profile around x
1
= 11 mm. When C
s
is in-
creased further to 400 GPa/m, a much smoother crack-like slip-rate profile is obtained.
Similarly, for the plastic calculations shown in Fig. 4.15(b), the train of pulses in the front
part of the slipping region is smoothed out and frictional slip propagates in a fully crack-
like mode with a non-slip region near the impact end.
Figure 4.16: Propagation speed of the frictional sliding front, V
slip
, versus time for the
elastic and plastic calculations with varying C
s
.
Fig. 4.16 shows that frictional slip initiates earlier and propagates slightly faster with
fewer oscillations for larger values of C
s
. Increased values of C
s
result in decreased val-
ues of
€
Δ˙ u
s
el
which means that
€
Δ˙ u
slip
€
Δ˙ u
s
as C
s
increases. Thus, the transition to a
smoother crack-like profile of
€
Δ˙ u
slip
with increasing C
s
is consistent with
€
Δ˙ u
s
tending to
remain smooth as noted in Section 4.3.2.
Although the mode of frictional slip is affected by varying C
s
, at least for the range of
values considered here, the partitioning of energy does not vary much.
0 5 10 15 20 25 30
1.5
2
2.5
3
3.5
4
Time (μs)
V
slip
(km/s)
c
l
C
s
= 100 GPa/m, elastic
C
s
= 100 GPa/m, plastic
C
s
= 200 GPa/m, elastic
C
s
= 200 GPa/m, plastic
C
s
= 400 GPa/m, elastic
C
s
= 400 GPa/m, plastic
95
In contrast to the effect of parameter C
s
, the dynamics of the frictional sliding are not
sensitive to variations in the parameter C
n
. We repeated the calculations for Case III us-
ing C
n
= 600 and 1200 GPa/m. The frictional slip propagation speed and the slip mode
are essentially unaffected by values of C
n
in this range. There are, however, some small
differences in the amplitude and phase of each pulse in the pulse train with different val-
ues of C
n
. For the elastic calculations, the pulses have slightly larger peak rates with
larger C
n
. For the plastic calculations, the pulse train near the rupture front is shorter with
larger values of C
n
. These small differences may be numerical artifacts since the adaptive
time step algorithm results in smaller times steps for larger values of C
n
.
4.3.4.2 Effect of impact velocity
To explore the effect of impact velocity, the elastic calculations and the plastic calcu-
lations with σ
0
= 10 MPa in Section 4.3.3 for Case III (where V
imp
= 20 m/s) are repeated
with V
imp
= 5 and 30 m/s.
In the elastic calculations, frictional slip propagates in a crack-like mode with V
imp
=
5 m/s (Fig. 4.17a), in contrast to propagating in a train-of-pulses mode with V
imp
= 20 m/s
(Fig. 4.11a) and 30 m/s (Fig. 4.18a). With V
imp
= 30 m/s the pulses are more closely
spaced and have higher peak frictional slip rates than in the corresponding calculation
shown in Fig. 11(a) with V
imp
= 20 m/s.
In the plastic calculation with V
imp
= 5 m/s as shown in Fig. 4.17(b), the profile of the
slip rate becomes very oscillatory and the frictional slip mode begins to develop into a
train of pulses. When V
imp
is increased to 30 m/s as shown in Fig. 4.18(b), the slip-rate
profile is much like that in Fig. 11(b) where V
imp
= 20 m/s.
96
(a) (b)
Figure 4.17: Along-interface profiles of
€
Δ˙ u
slip
and T
s
/T
n
at 10, 20 and 30 µs for (a) the
elastic calculation and (b) the plastic calculation with V
imp
= 5 m/s and all other parame-
ters the same as for Case III.
(a) (b)
Figure 4.18: Along-interface profiles of
€
Δ˙ u
slip
and T
s
/T
n
at 10, 20 and 28.5 µs for (a) the
elastic calculation and (b) the plastic calculation with V
imp
= 30 m/s and all other parame-
ters the same as for Case III.
Fig. 4.19 shows that frictional slip initiates earlier with increasing impact velocity for
both elastic and plastic calculations. In both cases, the slip propagation speeds with V
imp
=
20 and 30 m/s, approach essentially the same value, which is slightly lower than the cor-
responding calculations with V
imp
= 5 m/s. Thus, at least in the cases here, frictional slip
in a crack-like mode propagates at a slightly greater speed than propagation in a train-of-
pulses mode even though the impact velocity is smaller.
0 10 20 30 40 50 60 70
0
5
10
15
20
Δu
slip
(m/s)
x
1
(mm)
0
0.2
0.4
0.6
T
s
/
T
n
10 μs
20 μs
30 μs
0 10 20 30 40 50 60 70
0
20
40
60
Δu
slip
(m/s)
x
1
(mm)
0
0.2
0.4
0.6
T
s
/
T
n
10 μs
20 μs
30 μs
0 10 20 30 40 50 60 70
0
200
400
600
800
1000
Δu
slip
(m/s)
x
1
(mm)
0
0.2
0.4
0.6
T
s
/
T
n
10 μs
20 μs
28.5 μs
0 10 20 30 40 50 60 70
0
100
200
300
400
Δu
slip
(m/s)
0
0.2
0.4
0.6
T
s
/
T
n
10 μs
20 μs
28.5 μs
x
1
(mm)
97
Figure 4.19: Propagation speed of the frictional sliding front, V
slip
, versus time for the
elastic and plastic calculations with varying V
imp
.
(a) (b)
Figure 4.20: The effect of impact velocity V
imp
on the time variations of (a) the kinetic
energy W
kin
and (b) the frictional dissipation along the interface φ
slip
.
As seen in Fig. 4.20, there is an increasing difference between the predictions of the
elastic calculations and the plastic calculations for the energy partitioning with increasing
impact velocity: plastic dissipation increases, the kinetic energy of the plastic calculation
becomes increasingly less than that of the corresponding elastic calculation and the fric-
0 5 10 15 20 25 30 35
1
1.5
2
2.5
3
3.5
4
Time (μs)
V
imp
= 5 m/s, elastic
c
l
V
slip
(km/s)
V
imp
= 5 m/s, plastic
V
imp
= 20 m/s, elastic
V
imp
= 20 m/s, plastic
V
imp
= 30 m/s, elastic
V
imp
= 30 m/s, plastic
0
100
200
300
400
500
600
700
800
900
W
kin
(J/m)
V
imp
=
5
m/s, elastic
V
imp
=
5 m/s, plastic
V
imp
=
20
m/s, elastic
V
imp
=
20 m/s, plastic
V
imp
=
30
m/s, elastic
V
imp
=
30 m/s, plastic
0 5 10 15 20 25 30 35
Time (μs)
0 5 10 15 20 25 30 35
0
200
400
600
800
1000
1200
φ
slip
(J/m)
V
imp
=
5
m/s, elastic
V
imp
=
5 m/s, plastic
V
imp
=
20
m/s, elastic
V
imp
=
20 m/s, plastic
V
imp
=
30
m/s, elastic
V
imp
=
30 m/s, plastic
Time (μs)
98
tional dissipation of the plastic calculation becomes increasingly greater than that of the
corresponding elastic calculation.
4.4 Discussion
Two identical elastic-viscoplastic planar blocks having an initial compressive stress
were considered. Plane strain conditions were assumed and frictional sliding along the
interface was initiated by edge impact loading. Three sets of parameters were analyzed
that, for elastic material behavior, give rise to the same frictional slip modes as obtained
in the plane stress calculations of Coker et al. [2005], i.e., a crack-like mode, a pulse-like
mode and a train-of-pulses mode. In each case two sets of calculations were carried out,
an elastic calculation and a plastic calculation with a reference flow strength slightly
greater than the effective stress associated with the initial stress state. The effects of bulk
material plasticity on the frictional slip mode and on the energy partitioning were investi-
gated. The effect of variations in interface elastic properties and in impact velocity was
also explored.
With bulk plasticity, frictional slip initiates earlier but the slip propagation speeds are
not much affected. In the plastic calculations, the values of both T
s
and T
n
along the inter-
face are smaller than in the corresponding elastic calculations, and when the frictional
slip mode is crack-like, a zone of vanishing tractions along the interface develops from
the impact end of the interface.
99
(a)
(b)
Figure 4.21: Contours of effective plastic strain for (a) the plastic calculation in Case II
with σ
0
= 27 MPa at t = 37.8 µs and (b) the plastic calculation in Case III with σ
0
= 27
MPa at t = 30 µs.
In Case I (crack-like propagation) and in Case II (pulse-like propagation), the impact
velocity is relatively low, 2 m/s, and plastic deformation is rather limited (Fig. 4.21a).
The frictional slip mode and the energy partitioning are not much affected by bulk plastic
dissipation although the frictional slip rates are somewhat greater for plastic material be-
0.012
0.01
0.008
x
1
0 20 40 60 80
−40
−20
0
20
(mm)
x
2
(mm)
0.006
0.004
0.002
ε
ε
0.04
0.02
0.01
0.008
0.006
0.004
0.002
x
1
0 20 40 60 80
−40
−20
0
20
(mm)
x
2
(mm)
100
havior than for elastic material behavior. In Case III (train-of-pulses propagation), where
the impact velocity is 20 m/s and plastic deformation is more extensive (Fig. 4.21b), there
is a transition of the frictional slip mode from a train-of-pulses mode to a crack-like mode
in the calculation with bulk plasticity. Also, the more extensive plastic deformation in
this case results in an increase in frictional dissipation along the interface and a reduction
in kinetic energy as compared with the corresponding elastic calculation.
When the frictional slip rate
€
Δ˙ u
slip
has an abrupt change, the elastic slip rate
€
Δ˙ u
s
el
changes in the opposite direction with a shape and amplitude similar to
€
Δ˙ u
slip
; for exam-
ple when a pulse develops or when
€
Δ˙ u
slip
abruptly vanishes near the impact end. A con-
sequence of this is that the jump of the tangential velocity
€
Δ˙ u
s
maintains a relatively
smooth profile and is crack-like along the interface.
Calculations were carried out with the mesh size halved in each direction for a crack-
like case and for a pulse-like case with elastic bulk material behavior. For crack-like
propagation the results for the two meshes essentially coincided while for the pulse-like
case there was a small difference in the pulse profile along the interface. Calculations for
pulse-like frictional slip propagation with bulk material plasticity were also carried out
with initial time step values varying from 0.05 × min(ds)/c
l
to 0.001 × min(ds)/c
l
where
min(ds) is the minimum finite element edge length. The results agreed except for small
differences in the pulse amplitude in the early stages of slip.
The effects of varying the interface elastic parameters and the impact velocity were
explored and summarized in Fig. 4.22. Large values of C
s
tend to favor the crack-like
101
mode of sliding. In the plastic calculations, increased impact velocity gives rise to in-
creased bulk plastic deformation, increased frictional dissipation along the interface and a
reduction in kinetic energy. The effect of V
imp
on the frictional slip mode appears to be
more complicated. However, based on results for elastic calculations with C
s
= 200
GPa/m and plastic calculations with C
s
= 50 GPa/m, increasing V
imp
tends to favor a
crack-like mode or a train-of-pulses mode. Bulk material plasticity tends to smooth out
oscillations, which in turn tends to convert a train-of-pulses mode to a crack-like mode.
Figure 4.22: Phase diagrams showing the effects of impact velocity V
imp
and interface
stiffness parameter C
s
on the mode of frictional sliding with (a) elastic bulk material be-
havior and (b) plastic bulk material behavior. Shaded symbols denote the combinations of
parameters that produce different frictional sliding modes in elastic and plastic calcula-
tions. The initial compressive normal stress is Σ
0
= 10 MPa and the reference flow
strength for the plastic calculations is σ
0
= 10 MPa.
The relatively small effect of bulk material plasticity on the frictional slip mode and
on the propagation speed indicates that the interpretation of the Homalite on Homalite
2 5 10 20 30 50
50
100
200
400
V
imp
(m/s)
C
s
(GPa/m)
2 5 10 20 30 50
50
100
200
400
C
s
(GPa/m)
V
imp
(m/s)
(a) elastic calculations
(b) plastic calculations
Legend:
pulse-like mode a train of pulses or majority of Δu
slip
is a train of pulses over the range of calculation
crack-like mode but with notable oscillations in slip-rate profile
crack-like mode with relatively smooth slip-rate profile
102
sliding experiments in Coker et al. [2005] based on elastic material behavior likely re-
mains valid even if some bulk material plasticity did occur.
One of the main motivations for exploring the mechanics of fast frictional sliding
stems from earthquake dynamics. The loading conditions and bulk material constitutive
relation used in the calculations here are quite different from those appropriate for model-
ing earthquake ruptures. Whether or not the effects of plasticity seen in the circumstances
analyzed here carry over to the context of earthquake dynamics remains to be determined.
4.5 Appendix: a rate tangent method for updating the interface shear
traction
A rate tangent method [Povirk and Needleman, 1984] is used to update shear traction
rate
€
˙
T
s
as in Coker et al. [2005] and Shi et al. [2008]. The expressions used in the calcula-
tions are given here for completeness.
The frictional sliding rate,
€
Δ˙ u
slip
, in [t, t+Δt] is expressed as a linear combination of
its values at t and t+Δt
€
Δ˙ u
slip
= 1−γ
( )
Δ˙ u
slip
t
( )
+γΔ˙ u
slip
t +Δt
( ), (A4.1)
where 0 ≤ γ ≤ 1.
Expressing
€
Δ˙ u
slip
t +Δt
( ) as
€
Δ˙ u
slip
t +Δt ( ) =Δ˙ u
slip
t ( ) +Δt
∂Δ˙ u
slip
∂ T
s
sgn T
s
( )
˙
T
s
+
∂Δ˙ u
slip
∂θ
0
˙
θ
0
+
∂Δ˙ u
slip
∂θ
PC
˙
θ
PC
(A4.2)
and substituting into Eqs. (5.14) and (A.4.1) gives an expression of the form
€
˙
T
s
=− C
s
tan
Δ˙ u
s
−
˙
R
p [ ]
(A4.3)
103
with
€
C
s
tan
=
C
s
1+γΔt
∂Δ˙ u
slip
∂ T
s
C
s
(A4.4)
and
€
˙
R
p
= C
s
tan
sgn T
s
( ) Δ˙ u
slip
t ( ) +γΔt
∂Δ˙ u
slip
∂θ
0
˙
θ
0
+
∂Δ˙ u
slip
∂θ
PC
˙
θ
PC
. (A4.5)
The partial derivatives of
€
Δ˙ u
slip
in Eqs. (A4.4) and (A4.5) are obtained by direct dif-
ferentiation of Eq. (4.16)
€
∂Δ˙ u
slip
∂ T
s
=
rV
0
β
r−1
gθ
0
( )
θ
PC
, (A4.6)
€
∂Δ˙ u
slip
∂θ
0
=
−rV
0
β
r
gθ
0
( )
∂gθ
0
( )
∂θ
0
, (A4.7)
€
∂Δ˙ u
slip
∂θ
PC
=
−rV
0
β
r
θ
PC
, (A4.8)
where
€
∂g θ
0
( )
∂θ
0
=
µ
s
− µ
d
( )
P
θ
0
L
0
V
1
θ
0
p
exp −
L
0
V
1
θ
0
p
Q
1 r
+ µ
*
L
0
rV
0
θ
0
2
1
Q
1 r+1
(A4.9)
with
€
Q =
L
0
V
0
θ
0
+1, (A4.10)
€
µ
*
= µ
d
+ µ
s
− µ
d
( )exp −
L
0
V
1
θ
0
p
. (A4.11)
104
Chapter 5
Seismic Radiation From Tensile and Shear Point
Dislocations Between Similar and Dissimilar Sol-
ids
5.1 Introduction
Earthquakes are typically assumed to be caused by shear faulting [e.g., Aki and Rich-
ards, 2002]. While this assumption is made explicitly or implicitly in almost all earth-
quake-related studies, it is recognized that earthquake sources may also have tensile
components with differential motion in the fault-normal direction. Possible causes of
such fault-opening motion include volumetric changes in geothermal, volcanic and other
extensional regimes [e.g., Foulger and Long, 1984; Julian and Sipkin, 1985; Dreger et al.,
2000].
Recently there has been growing interest in the possibility that earthquake ruptures in
regular tectonic (not volcanic or geothermal) environments are associated with a small
amount of transient fault-opening motion. Such motion was observed in laboratory slid-
ing experiments between two foam rubber blocks [Brune et al., 1993] and earlier experi-
ments of Schallamach [1971] between rubbers and hard materials. Transient local fault
105
opening motion may be produced by collisions of geometrical asperities along rough sur-
faces [e.g., Lomnitz-Adler 1991], Schallamach-type waves, and dynamic changes of elas-
tic moduli in the source region [e.g., Knopoff and Randall, 1970; Ben-Zion and Ampu-
ero, 2009]. In addition, dynamic dilation near the tip of a rupture front propagating along
a bimaterial interface in the direction of slip on the compliant side of the fault can over-
come the initial compressive stress and produce local transient fault-opening [e.g.,
Weertman, 1980; Ben-Zion and Huang, 2002; Shi and Ben-Zion, 2006].
One obvious indication of fault-opening motion or tensile faulting may come from a
mismatch between fault-normal seismograms recorded at corresponding locations on the
opposite sides of a fault. However, such a comparison should take into account details of
the faulting geometry and the velocity structures on the opposite sides of the fault, and is
generally not easily achievable. Another set of signals that may reflect tensile faulting are
non-double-couple earthquake mechanisms [e.g., Julian et al., 1998; Miller et al., 1998].
However, non-double-couple mechanisms may also be produced by shear faulting in an
anisotropic medium [e.g., Kawasaki and Tanimoto, 1981] and a combination of different
double-couples [e.g., Frohlich et al., 1981; Kawakatsu, 1991; Bailey et al., 2009]. A third
type of signal, consistent with theoretical estimates of Haskell [1964] and Walter and
Brune [1993], are high-frequency S/P spectral amplitude ratios considerably smaller than
those expected for shear ruptures [e.g., Castro et al., 1991]. However, the high-frequency
spectral amplitudes are strongly affected by attenuation [e.g., Hanks, 1981] and the lack
of phase information does not allow detailed comparisons with various receiver positions.
106
Earthquake source studies usually also assume that the media on the opposite sides of
the fault have the same elastic properties. The assumed symmetry of properties across the
fault and the resulting symmetry of the radiation field from a shear dislocation lead to the
well-known problem of fault-plane ambiguity in focal mechanism studies using only di-
rect P-wave arrivals [e.g., Aki and Richards, 2002]. Many earthquakes, however, and es-
pecially moderate and large events occur on bimaterial interfaces that separate different
elastic solids. This is evident from tomographic images focusing on large faults [e.g.,
Eberhart-Phillips and Michael, 1998; Thurber et al., 2006] and analysis of head waves
that refract along bimaterial interfaces in the fault zone structure [e.g., Ben-Zion and Ma-
lin, 1991; Lewis et al., 2007]. The radiation pattern generated by sources on a bimaterial
interface should allow a distinction between the physical fault plane and the auxiliary
plane of the standard homogeneous models in theory [Ben-Zion, 1990].
To find additional potential seismic signatures of fault-opening motion and clarify the
radiation generated by sources on a bimaterial interface, we perform a theoretical study
involving analysis of synthetic seismograms generated by point sources with fault-
opening components between similar and dissimilar solids. Haskell and Thomson [1972]
and Thomson and Doherty [1977] calculated near-field seismograms from a finite propa-
gating tensile dislocation in an infinite homogeneous isotropic elastic solid. The wave-
form characteristics were affected by the finiteness and prescribed propagation speed of
the source. In the current study, however, we focus on detailed features of waveforms
produced by non-propagating localized shear and tensile dislocations. The results can be
107
considered as kernels for future calculations of seismic radiation from propagating finite
sources.
The synthetic waveform calculations indicate that near-fault seismograms for the
fault-parallel motion can be very useful for analysis of fault-opening motion, and that the
change of radiation (e.g., lose of symmetry and additional head wave phases) from
sources on bimaterial interfaces associated with moderate velocity contrast can be useful
for resolving the fault-plane ambiguity. The simulated seismograms provide guidelines
for observational studies aiming to derive high-resolution information on earthquake
source properties. The results can also be relevant for analysis of waves generated by ice-
quakes, ruptures in laboratory experiments, rock busts in mines, and seismic exploration
studies associated with explosions and fractures at the edges of boreholes.
5.2 Analysis
5.2.1 Radiation patterns from point sources in a homogeneous isotropic solid
Using the point-source approximation, we can derive the near-, intermediate-, and far-
field radiation patterns of body waves generated by a general moment tensor (Appendix
Section 5.4). The far-field patterns of the P- and S-waves from a shear dislocation in a
homogeneous solid can be expressed as [Aki and Richards, 2002]
€
A
S
FP
= sin2θ cosφ
( )
ˆ
r , (5.1a)
€
A
S
FS
= cos2θ cosφ ( )
ˆ
θ + −cosθ sinφ ( )
ˆ
φ (5.1b)
108
where
€
ˆ
r ,
€
ˆ
θ and
€
ˆ
φ are the unit basis vectors in spherical-polar coordinate system (Fig.
5.1). Similarly, the far-field radiation patterns of the P- and S-waves from a tensile dislo-
cation in a homogeneous solid (Appendix Section 5.4) are given by
€
A
T
FP
= λ µ + 2cos
2
θ
( )
ˆ r (5.2a)
€
A
T
FS
= −sin2θ ( )
ˆ
θ (5.2b)
where λ and µ are the Lamé parameters.
Figure 5.1: Cartesian and spherical-polar coordinate systems. The fault lies on the x
1
-x
2
plane with the x
3
axis as the fault normal. For a shear dislocation source the fault slips
along the x
1
-axis direction and for a tensile dislocation the fault displacement motion is
along the x
3
-axis direction.
Fig. 5.2 illustrates the far-field radiation patterns given by Eqs (5.1) and (5.2) for
shear and tensile dislocations. The assumed medium is a Poissonian whole space. In con-
trast to a shear dislocation, the radial and transverse components from a tensile disloca-
tion have no dependence on φ and therefore the radiation patterns are axial-symmetric. In
the case of a shear dislocation, there are two nodal planes for the P-wave, namely the
fault plane and the auxiliary plane, and no nodal plane for the S-wave. In the case of a
tensile dislocation, however, there is no nodal plane for the P-wave and the fault plane is
a nodal plane for the S-wave [e.g., Haskell, 1964]. Therefore, the problem of fault plane
x
1
x
2
x
3
ˆ r
θ
φ
θ ˆ
φ ˆ
0
109
ambiguity using only P-wave arrivals from a shear dislocation does not exist for the ten-
sile dislocation. Moreover, the radiation patterns from these two dislocation types reach
maximum amplitudes in different spatial directions. In the next two sections we examine
in more detail the radiation fields generated by shear and tensile dislocation sources
through calculations of synthetic waveforms.
Figure 5.2: Far-field radiation patterns of P-wave (top plots) and S-wave (bottom plots)
from (a) a shear dislocation source and (b) a tensile dislocation source in a Poissonian
solid (ν = 0.25). Fault-normal direction and slip directions are labeled and ‘N’ denotes
the neutral axis.
5.2.2 Synthetic seismograms for dislocation sources in a homogeneous solid
A numerical implementation of the 3D analytical solution of Ben-Zion [1990, 1999]
is employed to calculate synthetic seismograms generated by a general dislocation on an
interface separating similar (this section) or dissimilar (next section) elastic solids. The
solution is for a model of two joined half-spaces without a free surface. The model ge-
(a) (b)
110
ometry and the employed coordinate system are illustrated in Fig. 5.3(a). The fault lies on
the x-y plane and the dislocation source, located at the x-axis, can have fault-parallel
(shear in the x-direction) and/or fault-normal (tensile in the z-direction) components. A
ramp source-time function is used with a rise time of 0.01 second. In this section we set
the P-wave speed to be c
p
= 5196 m/s and the S-wave speed to be c
s
= 3000 m/s. The
source depth is fixed at x
s
= 1 km for all calculations done in this study. The receivers are
Figure 5.3: (a) A model consisting of two half-spaces with the x-y plane as the fault plane
and the z-axis in the fault-normal direction. The dislocation source (red star on the x-axis)
can have both fault-parallel (marked as 1) and fault-normal (marked as 2) components. x
s
,
z
r
and R are the source depth with respect to the receiver plane (fixed at 1 km), fault-
normal distance of the receiver and the hypocentral distance. (b) Various receiver con-
figurations on plane x = 0 employed in the calculations of synthetic seismograms: Con-
figuration ‘RCA’ with receivers along a straight line y = 0 which is normal to the fault
and through the epicenter; Configuration ‘RCB’ with receivers along two straight lines
parallel (and very close) to the fault trace with z = ±0.1 m; Configuration ‘RCC’ with re-
ceivers along a circle centered at the epicenter with radius r = 0.2 km; and Configuration
‘RCD’ with receivers along a straight line y = 1 km. The receivers in each configuration
are in pairs symmetric about the fault plane.
y
z
y
z
y
z
(b) (a)
RCB
RCC
RCA
z
x
y
1
2
Source
Receiver
Shear
Tensile
Fault plane
θ arccos =
()
z
r R
θ
R
x
s
z
r
RCD
y
z
111
located on the y-z plane in various configurations labeled in Fig. 5.3(b) as RCA, RCB,
RCC and RCD. The parameter θ defined in Fig. 5.3(a) is used in RCA and RCD as an
alternative indication for the fault-normal distance of the receivers.
Figure 5.4: Synthetic seismograms from a shear dislocation recorded with array RCD.
The red and blue traces are obtained from numerical implementation of the analytical so-
lution given by Ben-Zion [1990, 1999] whereas the black thin traces are obtained from
numerical implementation of the analytical solution given in Aki and Richards [2002].
As a part of code validation, we show in Fig. 5.4 comparisons between the synthetic
seismograms given by our code and numerical implementation of the analytical expres-
sions of body-wave radiations given in the Appendix (Section 5.4) for a shear dislocation.
Note that all the solutions here are for two half-spaces without free surface effect. The
receiver configuration RCD is employed. In this figure and later waveform results, blue
traces denote seismograms recorded at receivers with z < 0 (slower medium in the case of
dissimilar solids) and red traces denote seismograms recorded at receivers with z > 0
(faster medium in the case of dissimilar solids). We note that the amplitudes of the P and
S phases on every trace are well matched. Since the fault in this case is a nodal plane for
85
o
65
o
35
o
45
o
55
o
75
o
θ ( away from the fault)
Time (s)
Fault-normal Vertical
direct P
direct S
y
z
0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9
Fault-parallel
+
−
+
−
+
−
+
−
+
−
+
−
112
the P-wave and a maximum plane for the S-wave, the P-wave approaches zero and the S-
wave approaches maximum as the receiver gets closer to the fault, i.e., |θ| → 90°.
Figure 5.5: The radiation curves (amplitude profiles of P and S phases) obtained from
theoretical formula (solid lines) and measurements from synthetic seismograms (diamond
and circle symbols) from a shear dislocation in a homogeneous solid with (a) array RCA
where θ is defined in Fig. 5.3a, (b) array RCB and (c) array RCC where α is the angle
measured from +y axis clock-wise on the y-z plane.
Fig. 5.5 shows the amplitude profiles of P and S phases along different receiver con-
figurations from a shear dislocation. The solid lines are the theoretical radiation curves of
the P-wave (blue) or S-wave (red) amplitudes given by the sum of the intermediate- and
far-field radiation terms (Eqs A5.7-A5.10 in Appendix Section 5.4). The P and S ampli-
tudes in each panel are normalized by a common maximum factor. Since the P-wave and
S-wave are mixed in the near-field term (Eq. A5.6 in Appendix Section 5.4) and the con-
tribution from the near-field radiation is less than 5% of the total field at the receivers, the
near-field term is not included in plotting the solid lines. This might lead to some dis-
crepancies between the radiation curves and measurements from the whole-field synthetic
seismograms. For the receiver configuration RCA (Fig. 5.5a), the S-wave amplitude has a
local minimum at θ = ±45° and the P-wave amplitude reaches a maximum at θ ≈ ±55°.
For RCB (Fig. 5.5b), the receivers are close to the fault, which is a nodal plane for the P-
0 20 40 60 80 100 120 140 160 180
α
(clock-wise from +y axis)
0 10 20 30 40 50 60 70 80 90
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
Normalized Amplitude
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
Normalized Amplitude
θ
( away from the fault)
P theoretical
S theoretical
P synthetic
S synthetic
P theoretical
S theoretical
P synthetic
S synthetic
(a) (b) (c)
0
−2 −1 0 1 2 3
Along−strike distance (km)
Normalized Amplitude
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
0
P theoretical
S theoretical
P synthetic
S synthetic
−3
113
wave and maximum plane for the S-wave. The largest mismatch occurs for the S-wave
amplitudes near the ends of RCB (Fig. 5.5b).
Figure 5.6: The radiation curves (amplitude profiles of P and S phases) that are counter-
parts of plots in Fig. 5.5 from a tensile dislocation.
Fig. 5.6 gives corresponding amplitude profiles from a tensile dislocation. The P and
S amplitudes in Figs 5.6(a) and 5.6(c) are normalized as in Fig. 5.5 by a common maxi-
mum, while in Fig. 5.6(b) each of the P and S amplitude profiles is normalized with re-
spect to its own maximum due to the large difference in the amplitudes of these signals.
In this case, the largest discrepancies between the theoretical curves and results associ-
ated with the synthetic seismograms occur for the S-wave amplitudes at receivers close to
the fault (α = 0° and 180°) for the RCC array (Fig. 5.6c). The excellent general agreement
between the amplitude measurements from our synthetic seismograms and the theoretical
curves provides additional validation for the employed numerical code.
Fig. 5.7 shows the three-component synthetic seismograms from a shear dislocation
using receiver arrays RCA, RCB and RCC. For the receiver configuration RCA, the ver-
tical components are identically zero (Fig. 5.7a). For RCB both fault-parallel and vertical
components are several orders of magnitude smaller than the fault-normal components
(Fig. 5.7b). Fig. 5.8 shows the corresponding seismograms from a tensile dislocation. The
(a) (b) (c)
−2 −1 0 1 2 3
Along−strike distance (km)
Normalized Amplitude
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
P theoretical
S theoretical
P synthetic
S synthetic
−3
P theoretical
S theoretical
P synthetic
S synthetic
0 10 20 30 40 50 60 70 80 90
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
Normalized Amplitude
θ
( away from the fault)
0
P theoretical
S theoretical
P synthetic
S synthetic
0 20 40 60 80 100 120 140 160 180
α (clock-wise from +y axis)
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
Normalized Amplitude
114
vertical components for RCA (Fig. 5.8a) are identically zero and the fault-normal com-
ponents for RCB (Fig. 5.8b) are several orders of magnitude smaller than the other two
components.
In the case of a shear dislocation, the P and S arrivals at receivers across the fault
have the same polarities on the fault-normal components and opposite polarities on the
fault-parallel and vertical components (Fig. 5.7). In contrast, the situation is reversed in
the case of a tensile dislocation (Fig. 5.8). It is also interesting to note that at the receiver
configuration RCB the fault-parallel components are several orders larger than the fault-
normal components in the case of a tensile dislocation, while it is exactly the opposite in
the case of a shear dislocation.
115
Figure 5.7: Three-component synthetic seismograms from a shear dislocation in a homo-
geneous solid recorded with array (a) RCA (b) RCB and (c) RCC. The vertical compo-
nents in (a) are identically zero. Both fault-parallel and vertical components in (b) are
several orders of magnitudes smaller than the fault-normal components. The vertical
components in (c) are amplified by 10 times. The largest waveform amplitudes are la-
beled on the right side in unit of m.
0.2 0.3 0.4 0.5 0.6 0.7
85
o
65
o
35
o
45
o
55
o
75
o
θ ( away from the fault)
0.2 0.3 0.4 0.5 0.6 0.7 0.2 0.3 0.4 0.5 0.6 0.7
Time (s)
Fault-parallel Fault-normal Vertical
−1
0
1
0.5
−0.5
Time (s)
Fault-parallel Fault-normal Vertical
0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9
Along-strike distance (km)
Time (s)
Fault-parallel Fault-normal Vertical (x10)
0.1 0.2 0.3 0.4 0.5 0.1 0.2 0.3 0.4 0.5 0.1 0.2 0.3 0.4 0.5
180
o
135
o
90
o
45
o
0
o
α (clock-wise from +y axis)
(a)
(b)
(c)
Shear dislocation in a homogeneous solid
direct P direct S
direct P
direct S
direct P
direct S
2.68x10
−6
2.68x10
−6
2.57x10
−6
y
z
y
z
y
z
+
−
+
−
+
−
+
−
+
−
+
−
+
−
+
−
+
−
+
−
+
−
116
Figure 5.8: Three-component synthetic seismograms from a tensile dislocation source in
a homogeneous solid recorded with array (a) RCA (b) RCB and (c) RCC. The vertical
components in (a) and fault-normal components in (b) are identically zero. The largest
waveform amplitudes are labeled on the right side in unit of m.
Time (s)
Time (s)
(a)
(b)
(c)
Tensile dislocation in a homogeneous solid
direct P
direct S
direct P
direct S
direct S
0.2 0.3 0.4 0.5 0.6 0.7
85
o
65
o
35
o
45
o
55
o
75
o
θ ( away from the fault)
0.2 0.3 0.4 0.5 0.6 0.7 0.2 0.3 0.4 0.5 0.6 0.7
Time (s)
Fault-parallel Fault-normal Vertical
−1
0
1
0.5
−0.5
Fault-parallel Fault-normal Vertical
0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9
Along-strike distance (km)
Fault-parallel Fault-normal Vertical
0.1 0.2 0.3 0.4 0.5 0.1 0.2 0.3 0.4 0.5 0.1 0.2 0.3 0.4 0.5
180
o
135
o
90
o
45
o
0
o
α (clock-wise from +y axis)
2.04x10
−6
5.10x10
−7
1.06x10
−6
direct P
y
z
y
z
y
z
+
−
+
−
+
−
+
−
+
−
+
−
+
−
+
−
+
−
+
−
+
−
117
5.2.3 Synthetic seismograms for dislocation sources between dissimilar solids
To examine the effects of having dissimilar media across the fault on the seismic ra-
diations generated by shear and tensile sources, we use in this section a model consisting
of two dissimilar half-spaces. The P- and S-wave speeds of the faster medium (z > 0), c
p1
and c
s1
, are the same as in the previous section. The P- and S-wave speeds for the slower
medium (z < 0) , c
p2
and c
s2
, are set to be 10% smaller than those of the faster side: c
p2
=
c
p1
/1.1 = 5196/1.1 = 4723.6 m/s and c
s2
= c
s1
/1.1 = 3000/1.1 = 2727.3 m/s.
Figure 5.9: The radiation curves (amplitude profiles of P and S phases) measured from
synthetic seismograms for a shear dislocation between dissimilar solids (10% velocity
contrast) with array (a) RCA, (b) RCB and (c) RCC. Solid and dashed curves denote
measurements on the fast-medium side and slow-medium side, respectively.
Figure 5.10: The radiation curves (amplitude profiles of P and S phases) measured from
synthetic seismograms for a tensile dislocation between dissimilar solids (10% velocity
contrast) with array (a) RCA, (b) RCB and (c) RCC. Solid and dashed curves denote
measurements on the fast-medium side and slow-medium side, respectively.
P amp (fast)
S amp (fast)
P amp (slow)
S amp(slow)
P amp (fast)
S amp (fast)
P amp (slow)
S amp(slow)
P amp (fast)
S amp (fast)
P amp (slow)
S amp(slow)
30 40 50 60 70 80 90
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
Normalized Amplitude
θ
( away from the fault)
0
−2 −1 0 1 2 3
Along−strike distance (km)
Normalized Amplitude
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
0
−3 0 20 40 60 80 100 120 140 160 180
α (clock-wise from +y axis)
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
Normalized Amplitude
(a) (b) (c)
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
Normalized Amplitude
−2 −1 0 1 2 3
Along−strike distance (km)
Normalized Amplitude
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
0 20 40 60 80 100 120 140 160 180
α (clock-wise from +y axis)
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
Normalized Amplitude
(a) (b) (c)
30 40 50 60 70 80 90
θ
( away from the fault)
−3
P amp (fast)
S amp (fast)
P amp (slow)
S amp(slow)
P amp (fast)
S amp (fast)
P amp (slow)
S amp(slow)
P amp (fast)
S amp (fast)
P amp (slow)
S amp(slow)
118
Fig. 5.9 illustrates the amplitude profiles of P-waves (blue) and S-waves (red) for dif-
ferent receiver configurations from a shear dislocation along the bimaterial fault. Fig.
5.10 shows corresponding results from a tensile dislocation. The amplitude profiles on
the slower-medium and faster-medium sides are plotted with dashed and solid lines, re-
spectively. In both figures the P and S amplitudes in each case are normalized by a com-
mon maximum factor, except for Figs 5.9(b) and 5.10(b) where the P and S amplitudes
are normalized by the maximums of their own. In general, the wave amplitudes on the
slower-medium side are larger, as expected, than at corresponding positions on the faster-
medium side. Comparisons between Figs 5.5 and 5.6 with Figs 5.9 and 5.10 illustrate the
effects of the velocity contrast across the fault on the seismic radiation profiles. For ex-
ample, comparing Fig. 5.9(a) with Fig. 5.5(a) indicates that the S-wave maximum has
shifted from the fault plane (θ = ±90°) to around θ = ±75°, while its local minimum re-
mains around θ = ±45°. The most dramatic change occurs for the S-wave amplitude pro-
files from a shear dislocation at the array RCC (Fig. 5.9c): the minimum of the S-wave
amplitude profile in the homogeneous case (Fig. 5.5c) becomes the maximum and the
overall amplitude on the slower side of the fault is significantly larger than on the faster
side.
To get a more complete picture of the effects of dissimilar solids on the radiation pat-
terns, we plot in Figs 5.11 and 5.12 the P-wave and S-wave amplitudes generated at the
horizontal receiver plane x = 0 by a shear dislocation and a tensile dislocation, respec-
tively. A comparison between the left and right panels of Figs 5.11 and 5.12 highlights
the asymmetric radiation patterns across the fault in cases where the fault separates dis-
119
similar solids. Another distinctive feature in the radiation generated for two dissimilar
solids, seen clearly for both the P- and S-waves produced by both shear and tensile dislo-
Figure 5.11: The spatial distribution of P- and S-wave amplitudes on a receiver plane x =
0 radiated from a shear dislocation between (a) similar and (b) dissimilar solids. The am-
plitudes are normalized with the maximum value measured. The fault-normal distance
(along z-axis) and the along-strike distance (along y-axis) are both normalized by the
source depth x
s
. The solid line at x = 0 denotes the fault trace and the dot at (0, 0) denotes
the epicenter. The head waves are recorded as first arrivals in the region between the
dashed line and the solid line (fault).
cations, is the existence of lens-shaped areas with maximum wave amplitudes on the
slower-medium side within a fault-normal distance from the fault of z ≈ 0.5x
s
or θ ≈
63.4º. The lens-shaped areas with increased amplitude are generated by P or S head
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
P-wave Radiation
S-wave Radiation
−1.5 −1.0 −0.5 0 0.5 1.0 1.5
−1.5
−1.0
−0.5
0
0.5
1.0
1.5
Normalized along-strike distance (y/x
s
)
−1.5
−1.0
−0.5
0
0.5
1.0
1.5
Normalized along-strike distance (y/x
s
)
Normalized fault-normal distance (z/x
s
)
(a)
−1.5 −1.0 −0.5 0 0.5 1.0 1.5
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
P-wave Radiation
−1.5 −1.0 −0.5 0 0.5 1.0 1.5
−1.5
−1.0
−0.5
0
0.5
1.0
1.5
Normalized along-strike distance (y/x
s
)
S-wave Radiation
−1.5
−1.0
−0.5
0
0.5
1.0
1.5
Normalized along-strike distance (y/x
s
)
−1.5 −1.0 −0.5 0 0.5 1.0 1.5
Normalized fault-normal distance (z/x
s
)
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
(b)
120
waves contributing to the P or S direct arrivals when their duration times overlap. The
amplitude pattern within the lens-shaped area depends on the relations between the dif-
ferential arrival times of the head and direct waves and the rise time of the source-time
function. The differential arrival times of the head and direct waves decrease with in-
creasing fault-normal distance, up to a critical value beyond which head waves disappear
[e.g., Ben-Zion, 1990]. The pattern is also modified by the effect of geometrical spread-
ing.
Figure 5.12: The spatial distribution of amplitudes of P- and S-waves on the receiver
plane x = 0 radiated from a tensile dislocation between (a) similar and (b) dissimilar sol-
ids. The notations and symbols are the same in Fig. 5.11.
−1.5 −1.0 −0.5 0 0.5 1.0 1.5
−1.5
−1.0
−0.5
0
0.5
1.0
1.5
P-wave Radiation
0.5
0.6
0.7
0.8
0.9
1
S-wave Radiation
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
Normalized fault-normal distance (z/x
s
)
Normalized along-strike distance (y/x
s
)
(a) (b)
−1.5 −1.0 −0.5 0 0.5 1.0 1.5
−1.5
−1.0
−0.5
0
0.5
1.0
1.5
Normalized along-strike distance (y/x
s
)
P-wave Radiation
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
−1.5 −1.0 −0.5 0 0.5 1.0 1.5
−1.5
−1.0
−0.5
0
0.5
1.0
1.5
Normalized along-strike distance (y/x
s
)
S-wave Radiation
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
−1.5 −1.0 −0.5 0 0.5 1.0 1.5
−1.5
−1.0
−0.5
0
0.5
1.0
1.5
Normalized along-strike distance (y/x
s
)
Normalized fault-normal distance (z/x
s
)
121
Figure 5.13: Contours of the degree of expected interference between the head and direct
S waves based on analytical calculation of the differential travel-time between the differ-
ent phases, source rise time, and amplitude reduction due to geometrical spreading. The
notations and symbols are as in Fig. 5.11. See text for additional explanation.
Fig. 5.13 illustrates with analytical calculations the expected interference effect of the
head and direct waves for S waves radiation as a function of their differential travel-times
and the source rise time. For the region beyond the critical fault-normal distance desig-
nated in Fig. 5.13 by the dashed line, there is no head wave and hence no interference ef-
fect. For the region between the dash-dotted line and the fault in Fig. 5.13, the differential
travel-time is larger than 1.5 times the rise time, and hence there is little to no interfer-
ence effect. Between the dashed line and dash-dotted line, there is an overlap between the
head wave and the direct wave. The degree of interference between the two different
phases is represented by a cosine taper function with center and width equal to the rise
Normalized fault-normal distance (z/x
s
)
Normalized along-strike distance (y/x
s
)
122
time and 1.5 times the rise time, respectively, and amplitudes accounting for the geomet-
rical spreading of body waves.
Figs 5.14 and 5.15 show the three-component synthetic seismograms from shear and
tensile dislocations, respectively, with the receiver configurations RCA, RCB and RCC.
As expected, we observe P and S head waves in seismograms recorded on the slower side
of the fault within a critical distance z
c
=
€
r⋅tan cos
−1
c
2
c
1
( ) [ ]
, where r is the distance the
head wave travels along the bimaterial interface and c
1
and c
2
are the faster and slower P-
or S-wave speeds [Ben-Zion, 1990]. Due to the velocity difference of the solids, there is
an obvious time shift between the P- and S-wave phases at receivers on opposite sides of
the fault with equal hypocentral distances. Compared with the homogeneous case, seis-
mograms from a dislocation source between dissimilar solids have more complicated
waveform shapes, which masks somewhat the differences between the radiations gener-
ated by shear and tensile dislocation sources.
The largest waveform amplitudes at RCA and RCC for both shear and tensile disloca-
tions (Figs 5.14a, 5.14c, 5.15a, 5.15c) are ~1.5-2 times those of the homogenous case
(Figs 5.7a, 5.7c, 5.8a, 5.8c). On the other hand, for the RCB array the largest waveform
amplitudes (Figs 5.14b, 5.15b) remain at the same level or lower than those of the ho-
mogenous case (Figs 5.7b, 5.8b). It is also interesting to note that in the case of a shear
dislocation with array RCB (Fig. 5.14b), the amplitudes of the fault-parallel and vertical
components gained appreciably whereas the fault-normal components have slightly lower
amplitudes compared to the corresponding results in the homogeneous case (Fig. 5.7b).
In the case of a tensile dislocation and receiver array RCB, the fault-normal components
123
in the homogeneous case are several orders smaller than the other components (Fig.
5.8b). However, in the case of two dissimilar solids the amplitudes of the fault-normal
components with same source and receiver configuration (Fig. 5.15b) are comparable to
those of the other components.
124
Figure 5.14: Three-component synthetic seismograms from a shear dislocation between
two dissimilar solids (10% velocity contrast) recorded with array (a) RCA (b) RCB and
(c) RCC. P and S head waves are present on seismograms recorded within critical fault-
normal distances. The vertical components in (c) are amplified by 10 times. The largest
waveform amplitudes are labeled on the right side in unit of m.
Time (s)
Time (s)
(a)
(b)
(c)
Shear dislocation between two dissimilar solids
5.21x10
−6
direct P
direct S
direct P
direct S
direct P
direct S
Time (s)
0.2 0.3 0.4 0.5 0.6 0.7
85
o
65
o
35
o
45
o
55
o
75
o
θ ( away from the fault)
0.2 0.3 0.4 0.5 0.6 0.7 0.2 0.3 0.4 0.5 0.6 0.7
Fault-parallel Fault-normal Vertical
−1
0
1
0.5
−0.5
Fault-parallel Fault-normal Vertical
0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9
Along-strike distance (km)
Fault-parallel Fault-normal Vertical (x10)
0.1 0.2 0.3 0.4 0.5 0.1 0.2 0.3 0.4 0.5 0.1 0.2 0.3 0.4 0.5
180
o
135
o
90
o
45
o
0
o
α (clock-wise from +y axis)
1.30x10
−6
6.10x10
−6
direct P
head P
direct P
head P
y
z
y
z
y
z
+
−
+
−
+
−
+
−
+
−
+
−
+
−
+
−
+
−
+
−
+
−
125
Figure 5.15: Three-component synthetic seismograms from a tensile dislocation between
two dissimilar solids (10% velocity contrast) recorded with array (a) RCA (b) RCB and
(c) RCC. P and S head waves are present on seismograms recorded within critical fault-
normal distances. The vertical components in (c) are amplified by 10 times. The largest
waveform amplitudes are labeled on the right side in unit of m.
3.08x10
−6
Time (s)
Time (s)
(a)
(b)
(c)
Tensile dislocation between two dissimilar solids
direct P
direct S
head P
direct S
Time (s)
0.2 0.3 0.4 0.5 0.6 0.7
85
o
65
o
35
o
45
o
55
o
75
o
θ ( away from the fault)
0.2 0.3 0.4 0.5 0.6 0.7 0.2 0.3 0.4 0.5 0.6 0.7
Fault-parallel Fault-normal Vertical
−1
0
1
0.5
−0.5
Fault-parallel Fault-normal Vertical
0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9
Along-strike distance (km)
Fault-parallel Fault-normal Vertical (x10)
0.1 0.2 0.3 0.4 0.5 0.1 0.2 0.3 0.4 0.5 0.1 0.2 0.3 0.4 0.5
180
o
135
o
90
o
45
o
0
o
α (clock-wise from +y axis)
6.67x10
−7
2.95x10
−6
direct P
head P
direct S
head S
direct P
direct S
head S
head P
direct P
y
z
y
z
y
z
+
−
+
−
+
−
+
−
+
−
+
−
+
−
+
−
+
−
+
−
+
−
126
5.3 Discussion
Earthquake motions are dominated overall by tangential slip along faults. However, a
small superposed transient fault-opening motion can change dramatically the physics of
the processes involved and the energy partition between dissipation and radiation. As
mentioned in the introduction, many mechanisms including geometrical irregularities,
presence of fluids and contrasts of elastic properties across faults can produce a small
component of tensile motion. It is thus important to search for indicative signals of tran-
sient fault-opening motion that may be observable for faulting events dominated overall
by shear motion. Toward this end, we performed detailed calculations of synthetic seis-
mograms generated by shear and tensile sources between similar and dissimilar solids
using various source-receiver geometries.
The results indicate that careful examination of fault-parallel and fault-normal seis-
mograms at appropriate receiver configurations may provide evidence for the existence
of a small amount of fault-opening motion. Receivers located very close to the fault (e.g.,
array RCB) are expected to record essentially zero fault-parallel seismograms from a
shear dislocation (Fig. 5.7b). However, such receivers can have fault-parallel waveforms
with appreciable amplitudes (Fig. 5.8b) if the source has some tensile component. Near-
fault seismograms generated by a failure process consisting primarily of shear slip with a
small amount of superposed tensile motion will have considerably larger fault-normal
components than fault-parallel components. Differential fault-normal seismograms at
corresponding positions across the fault, and non-zero fault-parallel seismograms re-
corded very close to the fault, can indicate the existence of some tensile motion in the
127
source region. We also note that the P and S arrivals from a shear dislocation at receivers
on opposite sides of a fault in a homogeneous solid have the same polarities on the fault-
normal component, and opposite polarities on the fault-parallel and vertical components,
while for a tensile dislocation the situation is exactly reversed (Figs 5.7, 5.8). This may
not be useful for regular earthquakes for which shear faulting is dominant and likely to
mask the polarity signals. However, the discussed polarity information may be useful for
laboratory experiments and seismic exploration studies involving tensile cracking with
magnitude comparable to or larger than shearing.
The presence of velocity contrast across the fault makes noticeable differences in the
radiation patterns (Figs 5.5-5.6 and 5.9-5.12) and waveform characteristics (Figs 5.7-5.8
and 5.13-5.14) for both shear and tensile dislocation sources. In general, both the P and S
phases have larger amplitudes on the slower-medium side due to the lower elastic moduli
(Figs 5.11, 5.12). However, the velocity contrast across the fault masks the distinctions
between shear and tensile dislocations by changing the polarities of the first arrivals and
the amplitude ratios of the fault-parallel to fault-normal components of motion. Another
mechanism that can change the amplitude ratios of the fault-normal and fault-parallel
seismograms is supershear rupture near the vicinity of the receiver array [Aagaard and
Heaton, 2004]. However, the arrival times and additional phase characteristics can be
used easily to distinguish between supershear rupture, bimaterial interface and fault-
opening motion.
The discussed results are based on simple model calculations with localized non-
propagating sources between two similar or dissimilar uniform half-spaces without a free
128
surface. Observed seismograms are generally considerably more complex than the calcu-
lated waveforms, and include effects of finite propagating sources, structural heterogene-
ities off the fault, attenuations of high-frequency signals, etc.
Efforts to detect the calculated signals should employ high-resolution near-fault data
and analyses techniques such as stacking and de-convolution that amplify the target sig-
nals and reduce other effects. In some situations such as tensile ice-quakes, tensile crack-
ing in laboratory experiments, borehole studies with tensile fractures, and rock bursts in
mines with simple velocity structures between the sources and receivers, the results pro-
vide clear guidelines that can be used directly to identify tensile dislocation sources. For
natural earthquakes and other processes that are dominated by shear dislocations, the de-
tection of a tensile component of motion at the source remains a very challenging task. A
useful strategy should employ a combination of signals such as those discussed in the
present work and anomalously low S/P spectral ratio of the high-frequency radiated
waves [e.g., Haskell, 1964; Castrol et al., 1991; Walter and Brune, 1993], and isotropic
component of high-frequency radiation in regular tectonic environments [e.g., Spudich
and Chiou, 2008; Ford et al., 2009].
5.4 Appendix: derivation of radiation patterns from shear and tensile
point-dislocation sources
The body-wave radiation from a general moment tensor M corresponding to a point
source [Aki and Richards, 2002] can be expressed as
€
M
pq
*G
np,q
=u
n
N
t ( ) +u
n
IP
t ( ) +u
n
IS
t ( ) +u
n
FP
t ( ) +u
n
FS
t ( ), (A5.1)
129
where,
€
u
n
N
t
( )
=
15γ
n
γ
p
γ
q
−3γ
n
δ
pq
−3γ
p
δ
nq
−3γ
q
δ
np
4πρ
1
r
4
τM
pq
t−τ
( )
dτ
r α
r β
∫
, (A5.2)
€
u
n
IP
t
( )
=
6γ
n
γ
p
γ
q
−γ
n
δ
pq
−γ
p
δ
nq
−γ
q
δ
np
4πρα
2
1
r
2
M
pq
t−
r
α
, (A5.3a)
€
u
n
IS
t
( )
=−
6γ
n
γ
p
γ
q
−γ
n
δ
pq
−γ
p
δ
nq
−2γ
q
δ
np
4πρβ
2
1
r
2
M
pq
t−
r
β
, (A5.3b)
€
u
n
FP
t ( ) =
γ
n
γ
p
γ
q
4πρα
3
1
r
˙
M
pq
t−
r
α
, (A5.4a)
€
u
n
FS
t
( )
=−
γ
n
γ
p
−δ
np
( )
γ
q
4πρβ
3
1
r
˙
M
pq
t−
r
β
. (A5.4b)
Eqs (A5.2), (A5.3) and (A5.4) correspond to the near-field term, intermediate-field P
and S terms, and far-field P and S terms, respectively.
We assume isotropy condition for the medium in the source region and approximate
the moment tensor by taking average over the entire slip area Σ with
€
M
pq
= m
pq
⋅Σ, (A5.5)
where,
€
m
pq
=λu
k
ν
k
δ
pq
+ µ u
p
ν
q
+ u
q
ν
p ( )
and
€
u is the average slip vector. Here and here-
after Einstein summation convention is employed.
Substitution of Eq. (A5.5) into Eqs (A5.2)-(A5.4) yields
(1) Near-field term:
€
u
n
N
t
( )
=
µΣ
4πρr
4
A
N
⋅ u
p
t−τ
( )
⋅τdτ
r α
r β
∫
, (A5.6a)
where,
€
A
N
=30γ
n
γ
p
γ
q
ν
q
−6γ
n
ν
p
−6ν
n
γ
p
−6δ
np
γ
q
ν
q
. (A5.6b)
130
(2) Intermediate-field P term:
€
u
n
IP
t
( )
=
µΣ
4πρα
2
r
2
A
IP
⋅ u
p
t− r α
( ), (A5.7a)
where,
€
A
IP
=12γ
n
γ
p
γ
q
ν
q
+ λ µ−2 ( )γ
n
ν
p
−2ν
n
γ
p
−2δ
np
γ
q
ν
q
. (A5.7b)
(3) Intermediate-field S term:
€
u
n
IS
t
( )
=
µΣ
4πρβ
2
r
2
A
IS
⋅ u
p
t− r β
( ) , (A5.8a)
where,
€
A
IS
=−12γ
n
γ
p
γ
q
ν
q
−2γ
n
ν
p
−3ν
n
γ
p
−3δ
np
γ
q
ν
q ( )
. (A5.8b)
(4) Far-field P term:
€
u
n
FP
t
( )
=
µΣ
4πρα
3
r
A
FP
⋅
˙
u
p
t− r α
( ), (A5.9a)
where,
€
A
FP
=2γ
n
γ
p
γ
q
ν
q
+ λ µ ( )γ
n
ν
p
. (A5.9b)
(5) Far-field S term:
€
u
n
FS
t
( )
=
µΣ
4πρβ
3
r
A
FS
⋅
˙
u
p
t− r β
( ), (A5.10a)
where,
€
A
FS
=− 2γ
n
γ
p
γ
q
ν
q
−ν
n
γ
p
−δ
np
γ
q
ν
q ( )
. (A5.10b)
We then express the radiation pattern functions A
N
, A
IP
, A
IS
, A
FP
and A
FS
in the spheri-
cal-polar coordinate system with basis vectors
€
ˆ
r ,
€
ˆ
θ and
€
ˆ
φ (Fig. 5.1) for two end-member
dislocation sources: pure-shear and pure-tensile.
(1) Pure-shear dislocation source
In this case we have
€
ν
k
u
k
=0, (A5.11)
and the following correspondence
131
€
2γ
n
γ
q
ν
q
γ
p
u
p
↔ u sin2θ cosφ
ˆ
r { }
n
(A5.12a)
€
γ
n
ν
p
u
p
↔0 (A5.12b)
€
2ν
n
γ
p
u
p
↔ u sin2θ cosφ
ˆ
r − 2sin
2
θ cosφ
ˆ
θ
( ) { }
n
(A5.12c)
€
2γ
q
ν
q
u
n
↔ u sin2θ cosφ
ˆ
r + 2cos
2
θ cosφ
ˆ
θ − 2cosθ sinφ
ˆ
φ
( ) { }
n
(A5.12d)
where, {}
n
denotes the nth component of a vector.
Substitution of Eq. (A5.12) into Eqs (A5.6b), (A5.7b), (A5.8b), (A5.9b) and (A5.10b)
with the constraint of Eq. (A5.11) yields
€
A
S
N
= 9sin2θ cosφ ( ) ˆ r + −6cos2θ cosφ ( )
ˆ
θ + 6cosθ sinφ ( )
ˆ
φ (A5.13)
€
A
S
IP
= 4sin2θ cosφ ( ) ˆ r + −2cos2θ cosφ ( )
ˆ
θ + 2cosθ sinφ ( )
ˆ
φ (A5.14a)
€
A
S
IS
= −3sin2θ cosφ ( ) ˆ r + 3cos2θ cosφ ( )
ˆ
θ + −3cosθ sinφ ( )
ˆ
φ (A5.14b)
€
A
S
FP
= sin2θ cosφ
( )
ˆ
r (A5.15a)
€
A
S
FS
= cos2θ cosφ ( )
ˆ
θ + −cosθ sinφ ( )
ˆ
φ . (A5.15b)
Eqs (A5.13)-(A5.15) duplicate the classical results for shear dislocation given in Eq.
(4.33) of Aki and Richards [2002]. Below we derive corresponding radiation patterns
generated by a tensile dislocation.
(2) Pure-tensile dislocation source
In this case we have
€
ν
k
u
k
= u , (A5.16)
and the following correspondence
€
γ
n
γ
q
ν
q
γ
p
u
p
↔ u cos
2
θ
ˆ
r
{ }
n
(A5.17a)
132
€
γ
n
ν
p
u
p
↔ u
ˆ
r { }
n
(A5.17b)
€
2ν
n
γ
p
u
p
↔ u 2cos
2
θ
ˆ
r − sin2θ
ˆ
θ
( ) { }
n
(A5.17c)
€
2γ
q
ν
q
u
n
↔ u 2cos
2
θ
ˆ
r − sin2θ
ˆ
θ
( ) { }
n
(A5.17d)
where, {}
n
denotes the nth component of a vector.
Substitution of Eq. (A5.17) into Eqs (A5.6b), (A5.7b), (A5.8b), (A5.9b) and (A5.10b)
with the constraint of Eq. (A5.16) yields
€
A
T
N
= 18cos
2
θ− 6
( )
ˆ r + 6sin2θ ( )
ˆ
θ (A5.18)
€
A
T
IP
= λ µ− 2 + 8cos
2
θ
( )
ˆ r + 2sin2θ ( )
ˆ
θ (A5.19a)
€
A
T
IS
= 2− 6cos
2
θ
( )
ˆ r + −3sin2θ ( )
ˆ
θ (A5.19b)
€
A
T
FP
= λ µ + 2cos
2
θ
( )
ˆ r (A5.20a)
€
A
T
FS
= −sin2θ ( )
ˆ
θ (A5.20b)
The far-field radiation terms given by Eq. (A5.20) have similar forms to the radia-
tions given by Eq. (41) of Rice [1980].
In summary, the body-wave radiation patterns from a pure-shear dislocation source
are given by Eqs (A5.13)-(A5.15) and the corresponding patterns from a pure-tensile dis-
location are given by expressions (A5.18)-(A5.20). Graphic representations of these far-
field terms are given in Figs 5.2(a) and 5.2(b).
133
Bibliography
Aagaard, B. and Heaton, T., Near-source ground motions from simulations of sustained
intersonic and supersonic fault ruptures, Bull. Seismol. Soc. Am., 94, 2064-2078,
2004.
Abraham, F. F. and Gao, H. J., How fast can cracks propagate? Phys. Rev. Lett., 84,
3113–3116, 2000.
Achenbach, J. D. and Epstein, H. I., Dynamic interaction of a layer and a half-space. J.
Eng. Mech. Div., EM5, 27-42, 1967.
Adams, G. G., Self-excited oscillations of two elastic half-spaces sliding with constant
coefficient of friction, J. Appl. Mech., 62, 867-872, 1995.
Adams, G. G., Steady sliding of two elastic half-spaces with friction reduction due to in-
terface stick-slip, J. Appl. Mech., 65, 470-475, 1998.
Adams, G. G. An intersonic slip pulse at a frictional interface between dissimilar solids,
J. Appl. Mech., 68, 81-86, 2001.
Adda-Bedia, M. and Ben Amar, M., Self-sustained slip pulses of finite size between dis-
similar materials, J. Mech. Phys. Solids, 51, 1849-1861, 2003.
Aki, K. and Richards, P. G., Quantitative Seismology (2
nd
Edition), University Science
Books, 2002.
Ampuero, J. -P. and Ben-Zion, Y., Cracks, pulses and macroscopic asymmetry of dy-
namic rupture on a bimaterial interface with velocity-weakening friction, Geophys. J.
Int., 173, 674-692, 2008.
Ando, R. and Yamashita, T., Effects of mesoscopic-scale fault structure on dynamic
earthquake ruptures: Dynamic formation of geometrical complexity of earthquake
faults, J. Geophys. Res., 112, B09303, doi:10.1029/2006JB004612, 2007.
Andrews, D. J., A numerical study of tectonic stress release by underground explosions,
Seismol. Soc. Am. Bull., 63, 1375–1391, 1973.
Andrews, D. J., Rupture velocity of plane strain shear cracks, J. Geophys. Res., 81, 5679-
5687, 1976.
134
Andrews, D. J., Test of two methods for faulting in finite-difference calculation, Bull.
Seismol. Soc. Am., 89, 931-937, 1999.
Andrews, D. J., Rupture dynamics with energy loss outside the slip zone, J. Geophys.
Res., 110, B01307, 2005.
Andrews, D. J. and Ben-Zion, Y., Wrinkle-like slip pulse on a fault between different ma-
terials, J. Geophys. Res., 102, 553-571, 1997.
Anooshehpoor, A. and Brune, J. N., Wrinkle-like Weertman pulse at the interface be-
tween two blocks of foam rubber with different velocities, Geophys. Res. Lett., 23,
2025-2028, 1999.
Archuleta, R. J., A faulting model for the 1979 Imperial Valley earthquake, J. Geophys.
Res., 89, 4559-4585, 1984.
Bailey, I.W., Becker, T.W. and Ben-Zion, Y., Patterns of Co-seismic Strain Computed
from Southern California Focal Mechanisms, Geophys. J. Int., in press, 2009.
Beeler, N. M. and Tullis, T. E., Self-healing slip pulses in dynamic rupture models due to
velocity dependent strength, Bull. Seismol. Soc. Am., 86, 1130-1148, 1996.
Belytschko, T., Chiapetta, R. L. and Bartel, H. D., Efficient large scale nonlinear tran-
sient analysis by finite elements, Int. J. Numer. Methods Eng., 10, 579-596, 1976.
Ben-Zion, Y., The response of two half spaces to point dislocations at the material inter-
face, Geophys. J. Int., 101, 507-528, 1990.
Ben-Zion, Y., Corrigendum to “The response of two half spaces to point dislocations at
the material interface” by Ben-Zion (1990), Geophys. J. Int., 137, 580-582, 1999.
Ben-Zion, Y., Dynamic ruptures in recent models of earthquake faults, J. Mech. Phys.
Solids, 49, 2209-2244, 2001.
Ben-Zion, Y., Key Formulas in Earthquake Seismology, International handbook of
earthquake and engineering seismology Part B, pp. 1857-1875, Academic Press,
2003.
Ben-Zion, Y., A comment on “Material contrast does not predict earthquake rupture
propagation direction” by R. A. Harris and S. M. Day, Geophys. Res. Lett., 33,
L13310, 2006.
Ben-Zion, Y. and Ampuero, J. –P., Seismic radiation from regions sustaining material
damage, submitted to Geophys. J. Int., 2009.
135
Ben-Zion, Y. and Andrews, D. J., Properties and implications of dynamic rupture along a
material interface, Bull. Seismol. Soc. Am., 88, 1085-1094, 1998.
Ben-Zion, Y. and Huang, Y., Dynamic rupture on an interface between a compliant fault
zone layer and a stiffer surrounding solid, J. Geophys. Res., 107(B2), 2042, 2002.
Ben-Zion, Y., Katz, S. and Leary, P., Joint inversion of fault zone head waves and direct
P arrivals for crustal structure near major faults, J. Geophys. Res., 97, 1943-1951,
1992.
Ben-Zion, Y. and Malin, P., San Andreas Fault Zone Head Waves Near Parkfield, Cali-
fornia, Science, 251, 1592-1594, 1991.
Ben-Zion, Y., Peng, Z., Okaya, D., Seeber, L., Armbruster, J. G., Ozer, N. , Michael, A.
J., Baris, S. and Aktar, M.,
A shallow fault zone structure illuminated by trapped
waves in the Karadere-Duzce branch of the North Anatolian Fault, western Turkey,
Geophys. J. Int., 152, 699-717, 2003.
Ben-Zion, Y. and Rice, J. R., Slip patterns and earthquake populations along different
classes of faults in elastic solids, J. Geophys. Res., 100, 12959-12983, 1995.
Ben-Zion, Y. and Rice, J. R., Dynamic simulations of slip on a smooth fault in an elastic
solid, J. Geophys. Res., 102, 17771-17784, 1997.
Ben-Zion, Y. and Sammis, C. G., Characterization of fault zones, Pure Appl. Geophys.,
160, 677-715, 2003.
Ben-Zion, Y. and Shi, Z., Dynamic rupture on a material interface with spontaneous gen-
eration of plastic strain in the bulk, Earth Planet. Sci. Lett., 236, 486-496, 2005.
Beroza, G. C. and Mikumo, T., Short slip duration in dynamic rupture in the presence of
heterogeneous fault properties, J. Geophys. Res., 101, 22449-22460, 1996.
Bouchon, M. and Vallée, M., Observation of long supershear rupture during the Magni-
tude 8.1 Kunlunshan Earthquake, Science, 301, 824-826, 2003.
Brietzke, G. B. and Ben-Zion, Y., Examining tendencies of in-plane rupture to migrate to
material interfaces, Geophys. J. Int., 167, 807-819, 2006.
Broberg, K. B., Cracks and fracture, Academic Press, London, 1999.
Brune, J. N., Brown, S. and Johnson, P. A., Rupture mechanism and interface separation
in foam rubber models of earthquakes: A possible solution to the heat flow paradox
and the paradox of large overthrusts, Tectonophysics, 218, 59-67, 1993.
136
Burridge, R., Admissible speeds for plane-strain self-similar shear cracks with friction
but lacking cohesion, Geophys. J. R. Astron. Soc., 35, 439-455, 1973.
Caroli, C., Slip pulses at a sheared frictional viscoelastic/nondeformable interface, Phys.
Rew., E 62, 1729-1737, 2000.
Castro, R. R., Anderson, J. G. and Brune, J. N., Origin of high P/S spectral ratios from
Guerrero accelerograph array, Bull. Seismol. Soc. Am., 81, 2268-2288, 1991.
Cochard, A. and Madariaga, R., Complexity of seismicity due to highly rate dependent
friction, J. Geophys. Res., 101, 25321-25336, 1996.
Cochard, A. and Rice, J. R., Fault rupture between dissimilar materials: Ill-posedness,
regularization, and slip-pulse response, J. Geophys. Res., 105, 25891-25907, 2000.
Coker, D., Lykotrafitis, G., Needleman, A. and Rosakis, A. J., Frictional sliding modes
along an interface between identical elastic plates subject to shear impact loading, J.
Mech. Phys. Solids, 53, 884-922, 2005.
Dalguer, L. A., Irikura, K. and Riera, J. D., Simulation of tensile crack generation by
three-dimensional dynamic shear rupture propagation during an earthquake, J. Geo-
phys. Res., 108(B3), 2144, doi:10.1029/ 2001JB001738, 2003.
Das, S. and Aki, K., A numerical study of two-dimensional spontaneous rupture propaga-
tion, Geophys. J. R. Astr. Soc., 50, 643–668, 1977.
Das, S. and Kostrov, B. V., An investigation of the complexity of the earthquake source
time function using dynamic faulting models, J. Geophys. Res., 93, 8035-8050, 1988.
Day, S. M., Three-dimensional finite difference simulation of fault dynamics: Rectangu-
lar faults with fixed rupture velocity, Bull. Seismol. Soc. Am., 72, 705-727, 1982.
Day, S. M., Guang, Y. and Wald, D. J., Dynamic stress changes during earthquake rup-
ture, Bull. Seismol. Soc. Am., 88, 512-522, 1998.
Di Toro, G., Goldsby, D. L. and Tullis, T. E., Friction falls towards zero in quartz rock as
slip velocity approaches seismic rates, Nature, 427, 436-439, 2004.
Dieterich, J. H., The time-dependent friction and mechanics of stick-slip, Pure Appl.
Geophys., 116, 790-806, 1978.
Dieterich, J. H., Modeling of rock friction 1. Experimental results and constitutive equa-
tions, J. Geophys. Res., 84, 2161-2168, 1979.
137
Dieterich, J. H., Constitutive properties of faults with simulated gouge, Mechanical Be-
havior of Crustal Rocks, eds. N. L. Carter, M. Friedman, J. M. Logan, and D. W.
Stearns, Geophys. Monor. Ser., 24, AGU, Washington, D.C., pp. 103-120, 1981.
Dieterich, J. H., Earthquake nucleation on faults with rate- and state-dependent strength,
Tectonophysics, 211, 115-134, 1992.
Dieterich, J. H. and Kilgore, B., Implications of fault constitutive properties for earth-
quake prediction, Proc. Natl. Acad. Sci., USA 93, 3787-3794, 1996.
Dor, O., Ben-Zion, Y., Rockwell, T. K. and Brune, J. N., Pulverized rocks in the Mojave
section of the San Andreas Fault Zone, Earth Planet. Sci. Lett., 245, 642-654, 2006a.
Dor, O., Rockwell, T. K. and Ben-Zion, Y., Geologic observations of damage asymmetry
in the structure of the San Jacinto, San Andreas and Punchbowl faults in southern
California: A possible indicator for preferred rupture propagation direction, Pure
Appl. Geophys., 163, 301-349, 2006b.
Dor, O., Yildirim, C., Rockwell, T. K., Ben-Zion, Y., Emre, O., Sisk, M. and Duman, T.
Y., Geologic and geomorphologic asymmetry across the rupture zones of the 1943
and 1944 earthquakes on the North Anatolian Fault: possible signals for preferred
earthquake propagation direction, Geophys. J. Int., 173, 483-504, 2008.
Dreger, D. S., Tkalcic, H. and Johnston, M., Dilatational processes accompanying earth-
quakes in the Long Valley Caldera, Science, 288, 122-125, 2000.
Eberhart-Phillips, D. and Michael, A. J., Seismotectonics of the Loma Prieta, Californa,
region determined from three-dimensional V
p
, V
p
/V
s
, and seismicity, J. Geophys. Res.,
103, 21,099-21,120, 1998.
Festa, G. and Vilotte, J. -P., Influence of the rupture initiation on the intersonic transition:
Crack-like versus pulse-like modes, Geophy. Res. Lett., 33, L15320, 2006.
Ford, S. R., Dreger, D. S. and Walter, W. R., Identifying isotropic events using regional
moment tensor inversion, J. Geophys. Res., 114, B01306, 2009.
Foulger, G. and Long, R. E., Anomalous focal mechanisms: tensile crack formation on an
accreting plate boundary, Nature, 310, 43-45, 1984.
Freund, L. B., Dynamic Fracture Mechanics, Cambridge University Press, Cambridge,
1998.
Frohlich, C., Riedesel, M. A. and Apperson, K. D., Note concerning possible mechanisms
for non-double-couple earthquake sources, Geophys. Res. Lett., 16, 523-526, 1981.
138
Fuis, G. S., Clayton, R. W., Davis, P. M., Ryberg, T., Lutter, W. J., Okaya, D. A., Hauks-
son, E., Prodehl., C., Murphy, J. M., Benthien, M. L., Baher, S. A., Kohler, M. D.,
Thygesen, K., Simila, G. and Keller, G. R., Fault systems of the 1971 San Fernando
and 1994 Northridge earthquakes, southern California: Relocated aftershocks and
seismic images from LARSE II, Geology, 31, 171-174, 2003.
Goldsby, D. L. and Tullis, T. E., Low frictional strength of quartz rocks at subseismic slip
rate, Geophys. Res. Lett., 29, 1884, doi:10.1029/2002GL015240, 2002.
Hanks, T.C., The corner frequency shift, earthquake source models and Q, Bull. Seismol.
Soc. Am., 71, 597-612, 1981.
Harris, R. and Day, S., Effect of a low-velocity zone on a dynamic rupture, Bull. Seismol.
Soc. Am., 87, 1267-1280, 1997.
Haskell, N. A., Total energy and energy spectral density of elastic wave radiation from
propagating faults, Bull. Seismol. Soc. Am., 54, 1811-1841, 1964.
Haskell, N. A. and Thomson, K. C., Elastodynamic near-field of a finite, propagating ten-
sile fault, Bull. Seismol. Soc. Am., 63, 675-697, 1972.
Heaton, T. H., Evidence for and implications of self-healing pulses of slip in earthquake
rupture, Phys. Earth Planet. Inter., 64, 1-20, 1990.
Henry, C. and Das S., Aftershock zones of large shallow earthquakes: fault dimensions,
aftershock area expansion and scaling relations, Geophys. J. Int., 147(2), 272–293,
2001.
Husseini, M. I., Energy balance for motion along a fault, Geophys. J. R. Astro. Soc., 49,
699-714, 1977.
Jaeger, J. C. and Cook, N. G., Fundamentals of Rock Mechanics, Chapman and Hall, pp.
593, 1979.
Julian, B. R., Miller, A. D. and Foulger, G. R., Non-double couple earthquakes, 1. The-
ory, Rev. Geophys., 36, 525-549, 1998.
Julian, B. R. and Sipkin, S. A., Earthquake processes in the Long Valley caldera area,
California, J. Geophys. Res., 90, 11,155–11,169, 1985.
Kanamori, H., Anderson D. L. and Heaton, T. H., Frictional melting during the rupture of
the 1994 Bolivian Earthquake, Science, 279, 839-842, 1998.
Kawakatsu, H., Insignificant isotropic component in the moment tensor of deep earth-
quakes, Nature, 351, 50–53, 1991.
139
Kawasaki, I. and Tanimoto, T., Radiation patterns of body waves due to the seismic dis-
location occurring in an anisotropic source medium, Bull. Seismol. Soc. Am., 71, 37-
50, 1981.
Knopoff L. and Randall, M. J., The compensated linear vector dipole: a possible mecha-
nism for deep earthquakes, J. Geophys. Res., 75(26), 4957-4963, 1970.
Krieg, R. O. and Key, S. W., Transient shell response by numerical time integration, Int.
J. Numer. Methods. Eng., 7, 273-286, 1973.
Kubair, D. V., Geubelle, P. H. and Huang, Y. Y., Intersonic crack propagation in homo-
geneous media under shear-dominated loading: theoretical analysis, J. Mech. Phys.
Solids, 50, 1547-1564, 2002.
Lapusta, N., Modes of dynamic rupture on interfaces with nonlinear rate and state friction
laws, Proc. 11th Int. Conf. Fracture, Turin, Italy, 2005.
Lapusta, N., Rice, J. R., Ben-Zion, Y. and Zheng, G., Elastodynamic analysis for slow
tectonic loading with spontaneous rupture episodes on faults with rate- and state-
dependent friction, J. Geophys. Res., 105, 23765-23789, 2000.
Lewis, M. A., Ben-Zion, Y. and McGuire, J., Imaging the deep structure of the San An-
dreas Fault south of Hollister with joint analysis of fault-zone head and direct P arri-
vals, Geophys. J. Int., doi:10.1111/j.1365-246X.2006.03319.x, 2007.
Lewis, M. A., Peng, Z., Ben-Zion, Y. and Vernon, F. L., Shallow seismic trapping struc-
ture in the San Jacinto fault zone, Geophys. J. Int., 162, 867–881, 2005.
Linker, M. F. and Dieterich, J. H., Effects of variable normal stress on rock friction: Ob-
servations and constitutive equations, J. Geophys. Res., 97, 4923-4940, 1992.
Lockner, D. A. and Okubo, P. G., Measurements of frictional heating in granite, J. Geo-
phys. Res., 88, 4313-4320, 1983.
Lomnitz-Adler, J., Model for steady friction, J. Geophys. Res., 96, 6121-6131, 1991.
Lutter, W.J. et al., Upper Crustal Structure from the Santa Monica Mountains to the Si-
erra Nevada, Southern California: Tomographic Results from the Los Angeles Re-
gional Seismic Experiment, Phase II (LARSE II), Bull. Seismol. Soc. Am., 94, 619-
632, 2004.
Manighetti, I., Campillo, M., Sammis C., Mai P. M. and King, G., Evidence for self-
similar, triangular slip distributions on earthquakes: Implications for earthquake and
fault mechanics, J. Geophys. Res., 110, B05302, 2005.
140
Marone, C., Laboratory-derived friction laws and their application to seismic faulting,
Ann. Rev. Earth Planet. Sci., 26, 643-696, 1998.
McGarr, A., On relating apparent stress to the stress causing earthquake fault slip, J.
Geophys. Res. ,104, 3003– 3011, 1999.
McGuire, J. J. and Ben-Zion, Y., High-resolution imaging of the Bear Valley section of
the San Andreas Fault at seismogenic depths with fault-zone head waves and relo-
cated seismicity, Geophys. J. Int., 163, 152-164, 2005.
McGuire, J. J., Zhao, L. and Jordan, T. H., Predominance of unilateral rupture for a
global catalog of large earthquakes, Bull. Seismol. Soc. Am., 92, 3309-3317, 2002.
Michael, A. J. and Eberhart-Phillips, D., Relations among fault behavior, subsurface ge-
ology, and three-dimensional velocity models, Science, 253, 651-654, 1991.
Miller, A. D., Foulger, G. R. and Julian, B. R., Non-double couple earthquakes, 2. Obser-
vations, Rev. Geophys., 36, 551-568, 1998.
Needleman, A., A continuum model for void nucleation by inclusion debounding, J.
Appl. Mech., 54, 525-531, 1987.
Needleman, A., Numerical modeling of crack growth under dynamic loading conditions,
Comp. Mech., 19, 463-469, 1997.
Nielsen, S. B., Carlson, J. M. and Olsen, K. B., Influence of friction and fault geometry
on earthquake rupture, J. Geophys. Res., 105, 6069-6088, 2000.
Peng, Z., Ben-Zion, Y., Michael, A. J. and Zhu, L., Quantitative analysis of seismic
trapped waves in the rupture zone of the 1992 Landers, California earthquake: Evi-
dence for a shallow trapping structure, Geophys. J. Int., 155, 1021-1041, 2003.
Peirce, D., Shih, C. F. and Needleman, A., A tangent modulus method for rate dependent
solids, Comp. Struct., 18, 875-887, 1984.
Poliakov, A. N. R., Dmowska, R. and Rice, J. R., Dynamic shear rupture interactions
with fault bends and off-axis secondary faulting, J. Geophys. Res., 107(B11),doi:
10.1029/2001JB000572, 2002.
Povirk, G. L. and Needleman, A., Finite element simulations of fiber pull-out, J. Eng.
Mater. Tech., 115, 286-291, 1993.
Prakash, V., Frictional response of sliding interfaces subjected to time varying normal
pressures, J. of Tribology, 120, 97-102, 1998.
141
Prakash, V. and Clifton, R. J., Time-resolved dynamic friction measurements in pressure-
shear. Experimental Techniques in the Dynamics of Deformable Solids, Applied Me-
chanics Div., Vol. 165, American society of Mechanical Engineers, New York, pp.
33-48, 1993.
Ranjith, K. and Rice J. R., Slip dynamics at an interface between dissimilar materials, J.
Mech. Phys. Solids, 49, 341-361, 2001.
Ravi-Chandar, K., Lu, J., Yang, B. and Zhu, Z., Failure mode transitions in polymers un-
der high strain rate loading, Int. J. Fracture, 101, 33–72, 2000.
Rice, J. R., The mechanics of earthquake rupture, in Physics of the Earth’s Interior, eds.
A. M. Dziewonski and E. Boschi, Italian Physical Society / North Holland, Amster-
dam, pp. 555-649, 1980a.
Rice, J. R., Elastic wave emission from damage processes, J. Nondestr. Eval.,, 1(4), 215-
224, 1980b.
Rice, J. R., Spatio-temporal complexity of slip on a fault, J. Geophys. Res., 98, 9885-
9907, 1993.
Rice, J. R., New perspectives on crack and fault dynamics, in Mechanics for a new mil-
lennium, eds. H. Aref and J. W. Phillips, Kluer Academic Publications, pp. 1-23,
2001.
Rice, J. R., Lapusta N. and Ranjith, K., Rate and state dependent friction and the stability
of sliding between elastically deformable solids, J. Mech. Phys. Solids, 49, 1865-
1898, 2001.
Rice, J. R. and Ruina, A. L., Stability of frictional sliding, J. Appl. Mech., 50, 343-349,
1983.
Rice, J. R., Sammis, C. G. and Parsons, R., Off-fault secondary failure induced by a dy-
namic slip pulse, Bull. Seismol. Soc. Am., 95,109-134, 2005.
Rosakis, A. J., Samudrala, O. and Coker, D., Cracks faster than the shear wave speed,
Science, 284, 1337-1340, 1999.
Rosakis, A. J., Samudrala, O. and Coker, D., Intersonic shear crack growth along weak
planes, Mater. Res. Innov., 3, 236-243, 2000.
Rosakis, A. J., Intersonic shear cracks and fault ruptures, Adv. Phys., 51 (4), 1189-1257,
2002.
Rubin, A. and Gillard, D., Aftershock asymmetry/rupture directivity along central San
Andreas fault microearthquakes, J. Geophys. Res., 105, 19,095-19,109, 2000.
142
Ruina, A. L., Slip instability and state variable friction laws, J. Geophys. Res., 88,
10,359-10,370, 1983.
Samudrala, O., Huang, Y. and Rosakis, A. J., Subsonic and intersonic shear rupture of
weak planes with a velocity weakening cohesive zone, J. Geophys. Res., 107, B8 arti-
cle no. 2170, 2002.
Schallamach, A., How does rubber slide?, Wear, 17(4), 301-312, 1971.
Scholz, C. H., The Mechanics of Earthquakes and Faulting, Cambridge University Press,
pp. 471, 2002.
Scott, J. S., Masters, T. G and Vernon, F. L., 3-D velocity structure of the San Jacinto
fault zone near Anza, California - I. P waves, Geophys. J. Int., 119,611-626, 1994.
Shapiro, N. M., Campillo, M., Stehly, L. and Ritzwoller M. H., High-resolution surface-
wave tomography from ambient seismic noise, Science, 29, 1615-1617, 2005.
Shi, Z. and Ben-Zion, Y., Dynamic rupture on a bimaterial interface governed by slip-
weakening friction, Geophys. J. Int., 165, 469-484, 2006.
Shi, Z., Ben-Zion, Y. and Needleman, A., Properties of dynamic rupture and energy par-
tition in a solid with a frictional interface, J. Mech. Phys. Solids, 56, 5-24, 2008.
Spudich, P. and Chiou, B. S. J., Directivity in NGA earthquake ground motions: Analysis
using isochrone theory, Earthquake Spectra, 24, 279-298, 2008.
Templeton, E. L. and Rice, J. R., Off-fault plasticity and earthquake rupture dynamics: 1.
Dry materials or neglect of fluid pressure changes, J. Geophys. Res., 113, B09036,
doi:1029/2007JB005529, 2008.
Thomson, K. C. and Doherty, R., A correction to “Elastodynamic near-field of a finite,
propagating tensile fault” by Haskell and Thomson, Bull. Seismol. Soc. Am., 67,
1215-1217, 1977.
Thurber, C., Seismic tomography of the lithosphere with body waves, Pure Appl. Geo-
phys., 160, 717-737, 2003.
Thurber, C., Roecker, S., Zhang, H., Baher, S. and Ellsworth, W., Fine-scale structure of
the San Andreas fault and location of the SAFOD target earthquakes, Geophys. Res.
Lett., 31, L12S02, 2004.
Thurber, C., Zhang, H., Waldhauser, F., Hardebeck, J. Michael, A. and Eberhart-Phillips,
D., Three-dimensional compressional wavespeed model, earthquake relocations and
focal mechanisms for the Parkfield, California, Region, Bull. Seismol. Soc. Am., 96,
S38-S49, 2006.
143
Tsutsumi, A. and Shimamoto, T., High-velocity frictional characteristics of gabbro, Geo-
phys. Res. Lett., 24, 699-702, 1997.
Uenishi, K. and Rice, J. R., Universal nucleation length for slip-weakening rupture insta-
bility under nonuniform fault loading, J. Geophys. Res., 108, B1, 2042, 2003.
Walter, W. R. and Brune, J. N., Spectra of seismic radiation from a tensile crack, J. Geo-
phys. Res., 98, 4449-4459, 1993.
Weertman, J. J., Dislocations moving uniformly on the interface between isotropic media
of different elastic properties, J. Mech. Phys. Solids, 11, 197-204, 1963.
Weertman, J. J., Unstable slippage across a fault that separates elastic media of different
elastic constants, J. Geophys. Res., 85, 1455-1461, 1980.
Weertman, J. J., Subsonic type earthquake dislocation moving at approximately √2 shear
wave velocity on interface between half spaces of slightly different elastic constants,
Geophys. Res. Lett., 29(10), doi:10.1029/2001GL013916, 2002.
Xia, K., Rosakis, A. J. and Kanamori, H., Laboratory earthquake: the sub-Rayleigh-to-
supershear rupture transition, Science, 303, 1859-1861, 2004.
Xia, K., Rosakis, A. J., Kanamori, H. and Rice, J. R., Laboratory earthquakes along in-
homogeneous faults: directionality and supershear, Science, 308, 681-684, 2005.
Xu, X. -P., Needleman, A. and Abraham, F. F., Effect of inhomogeneities on dynamic
crack growth in an elastic solid, Modelling Simul. Mater. Sci. Eng., 5, 489-516, 1997.
Yamashita, T., Generation of microcracks by dynamic shear rupture and its effects on
rupture growth and elastic wave radiation, Geophys. J. Int., 143, 395-406, 2000.
Zhao, P., Peng, Z., Shi, Z., Lewis, M. A. and Ben-Zion, Y., Variations of the velocity
contrast and rupture properties of M6 Earthquakes along the Parkfield section of the
San Andreas Fault, submitted to Geophys. J. Int., 2009.
Zheng, G. and Rice, J. R., Conditions under which velocity-weakening friction allows a
self-healing versus a crack-like mode of rupture, Bull. Seismol. Soc. Am., 88, 1466-
1483, 1998.
Abstract (if available)
Abstract
This thesis examines dynamic ruptures along frictional interfaces and seismic radiation in models of earthquake faults separating similar and dissimilar solids with the goal of advancing the understanding of earthquake physics.
Linked assets
University of Southern California Dissertations and Theses
Conceptually similar
PDF
Interaction between dynamic ruptures and off-fault yielding characterized by different rheologies
PDF
Detailed properties of seismic waveforms and earthquake sources in Southern California
PDF
Symmetry properties, pulverized rocks and damage architecture as signatures of earthquake ruptures
PDF
Multi-scale imaging and monitoring of crustal and fault zone structures in southern California
PDF
Volumetric interactions between major ruptures and fault zones illuminated by small earthquake properties
PDF
Observations and modeling of dynamically triggered high frequency burst events
PDF
Multi-scale imaging of major fault zones in Southern California
PDF
Heterogeneity of earthquake stress drops, focal mechanisms and active fault zones
PDF
From laboratory friction to numerical models of fault dynamics
PDF
The impact of faulting and folding on earthquake cycles in collisional orogens
PDF
Analysis of waveform and catalog data of aftershocks for properties of earthquakes and faults
PDF
High-resolution imaging and monitoring of fault zones and shallow structures: case studies in southern California and on Mars
PDF
Applying automated techniques to large seismic datasets for systematic analyses of phases, source, and structure
PDF
Microseismicity, fault structure, & the seismic cycle: insights from laboratory stick-slip experiments
PDF
Determining fault zone structure and examining earthquake early warning signals using large datasets of seismograms
PDF
Characterizing fault behavior and earthquake surface expression on timescales of single events to multiple millennia
PDF
Seismic velocity structure and spatiotemporal patterns of microseismicity in active fault zones of Southern California
PDF
Integration and validation of deterministic earthquake simulations in probabilistic seismic hazard analysis
PDF
Towards an understanding of fault-system mechanics: from single earthquakes on isolated faults to millenial-scale collective plate-boundary fault-system behavior
PDF
Stress-strain characterization of complex seismic sources
Asset Metadata
Creator
Shi, Zheqiang
(author)
Core Title
Dynamic rupture processes and seismic radiation in models of earthquake faults separating similar and dissimilar solids
School
College of Letters, Arts and Sciences
Degree
Doctor of Philosophy
Degree Program
Geological Sciences
Degree Conferral Date
2009-08
Publication Date
08/04/2009
Defense Date
04/29/2009
Publisher
University of Southern California
(original),
University of Southern California. Libraries
(digital)
Tag
bimaterial interface,crack propagation,dynamic rupture,earthquake energy,Earthquakes,fault mechanics,friction,frictional sliding,numerical simulations,OAI-PMH Harvest,plasticity,rupture pulses,seismic radiation,supershear rupture
Language
English
Contributor
Electronically uploaded by the author
(provenance)
Advisor
Ben-Zion, Yehuda (
committee chair
), Becker, Thorsten W. (
committee member
), Jordan, Thomas H. (
committee member
), Nakano, Aiichiro (
committee member
), Sammis, Charles G. (
committee member
)
Creator Email
samzqshi@gmail.com,zheqians@usc.edu
Permanent Link (DOI)
https://doi.org/10.25549/usctheses-m2463
Unique identifier
UC1481831
Identifier
etd-Shi-3086 (filename),usctheses-m40 (legacy collection record id),usctheses-c127-179036 (legacy record id),usctheses-m2463 (legacy record id)
Legacy Identifier
etd-Shi-3086.pdf
Dmrecord
179036
Document Type
Dissertation
Rights
Shi, Zheqiang
Type
texts
Source
University of Southern California
(contributing entity),
University of Southern California Dissertations and Theses
(collection)
Repository Name
Libraries, University of Southern California
Repository Location
Los Angeles, California
Repository Email
cisadmin@lib.usc.edu
Tags
bimaterial interface
crack propagation
dynamic rupture
earthquake energy
fault mechanics
friction
frictional sliding
numerical simulations
plasticity
rupture pulses
seismic radiation
supershear rupture