of 26
Models of Stochastic, Spatially Varying Stress in the Crust Compatible
with Focal-Mechanism Data, and How Stress Inversions
Can Be Biased toward the Stress Rate
by Deborah Elaine Smith
*
and Thomas H. Heaton
Abstract
Evidence suggests that slip in earthquakes and the resultant stress changes
are spatially heterogeneous. If crustal stress from past earthquakes is spatially hetero-
geneous, then earthquake focal mechanisms should also be spatially variable. We
describe the statistical attributes of simulated earthquake catalogs, including hypocen-
ters and focal mechanisms, for a spatially
3D
, time-varying model of the crustal stress
tensor with stochastic spatial variations. It is assumed that temporal variations in stress
are spatially smooth and are primarily caused by plate tectonics. Spatial variations in
stress are assumed to be the result of past earthquakes and are independent of time for
periods between major earthquakes. It is further assumed that heterogeneous stress can
be modeled as a stochastic process that is specified by an autocorrelation function.
Synthetic catalogs of earthquake hypocenters and their associated focal mechanisms
are produced by identifying the locations and times at which the second deviatoric stress
invariant exceeds a specified limit. The model produces a seismicity catalog that is
spatially biased. The only points in the grid that exceed the failure stress are thosewhere
the heterogeneous stress is approximately aligned with the stress rate. This bias results
in a focal-mechanism catalog that appears less heterogeneous than the underlying stress
orientations. Comparison of synthetic focal-mechanism catalogs with catalogs of real
earthquakes suggests that stress in the crust is heterogeneous. Stochastic parameters are
estimated which generate distance dependent spatial variations in focal mechanisms
similar to those reported by
Hardebeck (2006)
for southern California.
Introduction
Statistical Paradigm for Stress Heterogeneity
Seismically active regions in the Earth have a great deal
of complexity that can lead to spatially heterogeneous stress.
Specifically, in many of these seismic regions, the crust is
highly fractured by faults of varying lengths. Then the faults
within these networks can have ruptures with spatially vari-
able slip. Given this complex spatial geometry and slipping
histories, integrated through time for all events, it follows
that the stress distribution should also be highly complex
or spatially heterogeneous.
A natural question therefore arises: what is the best char-
acterization of heterogeneous stress in the crust for these
regions? Given the complex slip history of the Earth
s crust
at widely varying scales, it is not reasonable to determinis-
tically track stress changes from every dynamic event over
every length scale. At the same time, the development of
a heterogeneous stress model in
3D
that encompasses all
length scales would be a significant improvement over com-
monly used homogeneous stress models. Therefore, as a
first-order approximation, we construct the following simple
3D
spatially heterogeneous stress model. We assume that the
spatially heterogeneous stress is independent of time for the
interseismic period (between major events) and that it can be
described with two parameters that characterize fractal-like
statistics for each component of the
3D
stress tensor. This
fractal-like spatial variation of stress represents the stress
pattern left over from the combined action of previous earth-
quakes. We also assume that during the interseismic period,
stress from tectonic forces grows steadily in time and
smoothly in space. This steady stress change is in turn super-
imposed on the rough fabric of stress heterogeneity. This
simple model captures spatial variations of stress at all length
scales, has two parameters that allow us to tune the fractal-
like statistics of the stress to better match data, and has load-
ing that can bring points to failure given a fracture criterion.
Another question we investigate with this simple
model is, What effect can stress heterogeneity have on stress
*Now at Department of Terrestrial Magnetism, Carnegie Institution of
Washington, 5241 Broad Branch Road NW, Washington, D.C. 20015-1305.
1396
Bulletin of the Seismological Society of America, Vol. 101, No. 3, pp. 1396
1421, June 2011, doi: 10.1785/0120100058
inversion estimates of background stress (the spatial mean
stress within a region)? In particular, we show with our
model that if stress is spatially variable, the hypocenters of
events tend to occur only at those locations where the spa-
tially heterogeneous shear stress is substantially increasing
with time. In other words, ruptures are favored at points
where the heterogeneous stress and stressing rate are aligned,
and ruptures are suppressed when the heterogeneous stress
and increasing stress are misaligned. This means that earth-
quake focal mechanisms can represent a biased or nonuni-
form sampling of the spatially heterogeneous stress field
because some stress orientations are favored for failure over
others depending upon their relative orientation to the stres-
sing rate. This biased sampling of the stress field by focal
mechanisms in turn produces biased stress inversion results.
Instead of reproducing the spatial mean stress orientations
within a region, the stress inversions now generate principal
stress orientations biased toward the stress rate (which may
not be aligned with the spatial mean).
Hence, our goal in this paper is to develop stochastic
models of
3D
deviatoric stress heterogeneity that: (1) have
first-order realism; that is, the models can generate focal
mechanisms statistics matching those seen in the real Earth
and (2) help parameterize the potential stress inversion bias
toward the stressing rate orientation. To do this, we employ
3D
numerical grids of dimensions
201
×
201
×
201
points
and
2D
numerical grids of dimensions
3375
×
3375
points,
with stress defined at each point and as a function of time.
We define the stress in these grids with only two statis-
tical parameters,
α
, which is related to the fractal dimension
or spatial roughness of stress, and heterogeneity ratio (
HR
),
which describes the ratio of spatial heterogeneity to the spa-
tial mean. The stress is defined at each point as a function of
time in
3D
numerical grids of dimensions
201
×
201
×
201
points and in
2D
numerical grids of dimensions
3375
×
3375
points. We generate synthetic earthquakes and their asso-
ciated focal mechanisms from these numerical models by
applying an appropriate failure criterion at each point. We
will show that the different values of
α
and
HR
have distinct
effects on the synthetic seismicity catalogs. By comparing
the orientation and clustering statistics of our synthetic
data sets with a catalog of focal mechanisms in southern
California, we produce an estimate of our statistical stress
parameters
α
and
HR
. This comparison also matches a 60 km
estimate to the outer scale of our grids; hence, we have a
maximum resolution of about 300 to 600 m for the
3D
grids
and 20 to 40 m in
2D
. We then use these synthetic data sets,
with a focus on the preferred
α
and
HR
, to estimate stress
inversion orientation biasing for southern California.
Our model to create heterogeneous focal-mechanism
orientations contrasts with the stress model used by many
previous studies of focal-mechanism statistics. In these pre-
vious studies, stress is assumed to be approximately homo-
geneous, and failure is assumed to occur on randomly
oriented planes of weakness (preexisting faults), which re-
quires variable frictional strength. Consequently, the studies
invert for a spatially homogenous stress that best aligns with
potential slip vectors from focal-mechanism catalogs (
Carey
and Brunier, 1974
;
Angelier, 1975
;
Etchecopar
et al.
, 1981
;
Angelier, 1984
;
Gephart and Forsyth, 1984
;
Michael, 1984
,
1987
;
Mercier and Carey-Gailhardis, 1989
;
Gephart, 1990
).
While these homogeneous stress models can apparently
explain many features of focal-mechanism characteristics,
we will demonstrate that it is also possible to explain these
statistics with our model, which is quite different. In our
model, the stress has large amplitude variations over short
length scales (which is consistent with spatial variations in
fault slip), homogeneous temporally varying stress due to
plate tectonic loading, and uniform frictional strength.
It seems clear that these two classes of models (homo-
geneous stress and heterogeneous strength versus heteroge-
neous stress and homogeneous strength) represent end
members, and that the real Earth likely has both heteroge-
neous stress and heterogeneous strength (
Rivera and Kana-
mori, 2002
). While each of these models allows researchers
to infer Earth processes from focal-mechanism catalogs, we
will show that the model of heterogeneous stress and homo-
geneous strength provides very different information than
can be obtained from traditional stress inversion studies. It
is important to explore the heterogeneous stress model and
its consequences.
While there are previous studies that attempted to exam-
inetheeffectofheterogeneous stress(
Michael,1991
;
Lu
etal.
,
1997
) on stress inversion results, they do not resolve the issues
raised in this paper for highly heterogeneous stress. For
example, the spatially heterogeneous component of stress
in our models, with a best fit to data, is distinctly larger than
the spatially uniform component of stress. We also show that
to determine the effect of stress heterogeneity, it is critical that
the stress rate is accounted for concurrently with stress hetero-
geneity. It is the interaction between the stress heterogeneity
and the stress perturbation term (in this case the stressing rate)
that produces a biased sampling of which points fail and are
included in the synthetic focal-mechanism catalogs.
Evidence for Spatially Heterogeneous Stress
Observations of spatially varying slip along fault zones
and in earthquakes suggest that both slip and stress are very
spatially heterogeneous and possibly fractal in nature
(
Andrews, 1980
,
1981
;
Herrero and Bernard, 1994
;
Manigh-
etti
et al.
, 2001
;
Mai and Beroza, 2002
;
Ben-Zion and Sam-
mis, 2003
;
Lavallee and Archuleta, 2003
;
Manighetti
et al.
,
2005
). For example,
McGill and Rubin (1999)
observed a
1 m change in slip over distances of approximately 1 km in
the Landers earthquake, which indicates a
10

3
strain change.
This implies possibly a 100 MPa stress change averaged over
the distance of 1 km. Similar strain changes can be seen in the
slip inversion of seismic and geodetic data from the Landers
earthquake (
Wald and Heaton, 1994
).
Another example of highly variable, heterogeneous slip
over short wavelengths comes from
Manighetti
et al.
(2001)
.
Models of Stochastic, Spatially Varying Stress in the Crust Compatible with Focal-Mechanism Data
1397
Using altimetry data in the Afar depression, East African rift,
they show heterogeneous cumulative slip as a function of
distance, with short wavelength strains of the order
5
×
10

2
.
While nonelastic processes may come into play at such large
shear strains, these observations of heterogeneous slip
demonstrate certain features. Heterogeneous slip patterns
exist not just for individual earthquake slip histories, but they
persist for the entire cumulative slip history of fault zones,
indicating that slip heterogeneity is a stable feature. In addi-
tion, the cumulative slip shows possibly self-similar, fractal
patterns (
Manighetti
et al.
, 2001
;
Manighetti
et al.
, 2005
).
Borehole studies, which measure the orientation of
maximum horizontal compressive stress,
S
H
, directly from
borehole breakouts, also indicate that stress can be quite
heterogeneous.
Wilde and Stock (1997)
reported on multiple
boreholes with different orientations that had been drilled at
approximately the same locations. They analyzed this data to
constrain the relative magnitudes of the principal stresses.
Boreholes drilled within close proximity of each other (less
than 1 km) show greatly varying
S
H
orientations (tens of
degrees), which may be indicative of heterogeneous stress
(see Fig.
1
). The orientations of breakouts in the Cajon Pass
borehole (
Barton and Zoback, 1994
) also show significant
heterogeneity for an individual borehole near an active fault.
Figure
2
shows some of this orientation data. Later we will
compare it side by side with synthetically generated orienta-
tion data from our stochastic stress models in
The Size of
Stress at Different Scales
.
Liu-Zeng
et al.
(2005)
have also shown that the assump-
tion of short wavelength heterogeneous fractal slip can repro-
duce distributions of earthquakes having slip versus length
ratios similar to real earthquakes and realistic Gutenberg
Richter frequency magnitude statistics. Using simple
stochastic models, they showed that spatially connected slip
can produce average stress drops (a constant times average
slip divided by rupture length) similar to real data.
Perhaps the most interesting piece of data comes from
Zoback and Beroza (1993)
. They studied the orientations of
aftershock planes from the Loma Prieta earthquake and
plotted their distributions as a function of strike and dip.
They found aftershocks that had both right-lateral and left-
lateral orientations on similar fault planes as well as normal
and reverse orientations. Given that this seismicity is in the
immediate vicinity of the right-lateral San Andreas fault, the
existence of left-lateral aftershocks on fault planes parallel to
the San Andreas fault presents a curious problem. Zoback
and Beroza proposed that the principal compressive stress
direction was almost normal to the San Andreas fault and
that the aftershocks occurred on extremely weak faults of
different orientations surrounding the mainshock zone.
However, if one allows for a paradigm of spatially heteroge-
neous stress in three dimensions, as presented in this paper,
the left-lateral orientations naturally occur.
The model we present predicts that such opposite
mechanisms should primarily be observed for aftershocks.
Figure
3
shows our hypothesis for what a
1D
cross section
of shear stress in southern California might look like. While
most of the points have positive shear stress on the
σ
0
12
plane,
a fair percentage have negative shear stress on the
σ
0
12
plane.
Heterogeneity similar to this could explain why Zoback and
Beroza observed left-lateral aftershocks after the Loma Prieta
earthquake. In particular, the large local stress change to the
Figure
1.
Wilde and Stock (1997)
plotted inferred maximum horizontal compressive stress,
S
H
, orientations from borehole breakouts in
southern California. There are a variety of orientations for borehole breakouts from the same borehole or from boreholes spatially close to one
another. This suggests short-wavelength spatial stress heterogeneity. In this modified plot, we use circles to point out a few of the locations
studied by Wilde and Stock that show evidence for
S
H
orientation heterogeneity. The color version of this figure is available only in the
electronic edition.
1398
D. E. Smith and T. H. Heaton
system from the mainshock, combined with stress heteroge-
neity in the left-lateral direction, would create left-lateral
aftershocks.
Indeed, others have also suggested that spatially hetero-
geneous stress could explain the Loma Prieta observations.
Gephart (1997)
studied the first six weeks of aftershock data
in six subregions and found substantial variation in the
inferred stress directions from earthquake focal mechanisms.
Michael
et al.
(1990)
suggested that the Loma Prieta after-
shocks are indicative of highly heterogeneous stress, and this
heterogeneous stress (with implied localized stress drops of
100 MPa or more) could be created by irregular slip in
the rupture. Two papers in particular (
Eberhart-Phillips
and Michael, 1998
;
Michael and Eberhart-Phillips, 2000
)
demonstrate how the Loma Prieta aftershock data are more
consistent with a postseismic heterogeneous stress field than
the Zoback and Beroza conclusion of a weak fault with near
perpendicular maximum horizontal compressive stress.
Stress Model to Be Used
We develop the following fractal-like model of crustal
stresses in space and time. We assume that the fractal-like
statistics are approximately time independent; hence, this
description is intended only for interseismic times, the peri-
ods of time between large earthquake sequences. In our
numerical models, we construct
3D
and
2D
spatial grids
with the deviatoric stress defined at each grid point, where
σ
0
ij

σ
ij

δ
ij
σ
kk
=
3
denotes deviatoric stress. Deviatoric
stress is sufficient for our simulations because the fracture
criterion we use (see
Fracture Criterion Used to Create
Figure
2.
On the left is borehole breakout data from the Cajon
Pass in a particularly heterogeneous section (
Barton and Zoback,
1994
). In this figure, modified from
Barton and Zoback (1994)
, plus
signs are used to plot the maximum horizontal compressive stress,
S
H
, azimuth as a function of depth and triangles are used to plot
their model. We added the arrows to show that
S
H
orientations
can easily vary by 90°. The color version of this figure is available
only in the electronic edition.
0
20
40
60
80
100
−100
−50
0
50
100
150
200
Length (km)
Shear Stress (MPa)
Sample 1D Synthetic Shear Stress Profile
Figure
3.
Mean shear stress for a synthetic
1D
profile, generated
by the model presented in this paper. Even though the mean shear
stress over the entire length is quite positive, the spatial variation
creates pockets of shear stress with the opposite sign. This type
of model could explain why Loma Prieta had some left-lateral
aftershocks on a predominantly right-lateral fault system. The color
version of this figure is available only in the electronic edition.
Models of Stochastic, Spatially Varying Stress in the Crust Compatible with Focal-Mechanism Data
1399
Synthetic Focal-Mechanism Catalogs
) to generate synthetic
focal mechanisms is purely a function of deviatoric stress.
The total deviatoric stress is constructed as follows:
σ
0

x
;t

σ
0
B

_
σ
0
T

t

t
0

σ
0
H

x

:
(1)
The term
σ
0
B
is the background stress, which is the spatially
and temporally averaged deviatoric stress tensor in the
region of interest. The intended solution of stress inver-
sions is
σ
0
B

_
σ
0
T

t

t
0

minus any overall stress magni-
tude information. For the values used in this paper,
σ
0
B

_
σ
0
T

t

t
0

σ
0
B
; therefore, we are interested in seeing how
well stress inversions can resolve the three principal orienta-
tions of
σ
0
B
and the stress ratio of
σ
0
B
,
R

σ
0
B
2

σ
0
B
3

=

σ
0
B
1

σ
0
B
3

(
Rivera and Kanamori, 2002
).
The term
_
σ
0
T

t

t
0

is the temporally varying deviatoric
stress due to plate tectonics that brings points within our
3D
grid to failure as point earthquakes.
t
0
is the time of the last
large earthquake where
t
0
is arbitrarily chosen to be the start
time of any of our simulations. We assume that temporal
variations in stress are primarily caused by forces that are
applied at a distance and that the temporally varying stress
is approximately spatially homogeneous. While it is clear
that temporal variations in stress do change with location,
the observed spatial variations are small compared with
spatial heterogeneities that arise from heterogeneous slip in
past earthquakes.
This term is assumed to grow linearly with time for our
simulation time windows, where we choose a magnitude of
1
MPa
=
century for
_
σ
0
T
. We derive this approximate magni-
tude from estimates of the strain rate in southern California
and especially near the San Andreas fault system of the order
0
:
2
μ
strain
=
yr (
Johnson
et al.
, 1994
;
Shen
et al.
, 1996
),
coupled with shear modulus estimates in the range of 40 GPa
(
Turcotte and Schubert, 1982
).
The time windows are set by the time it takes for the first
2000 points within our grid to fail; therefore, for a
_
σ
0
T
of
magnitude
1
MPa
=
century, the first 2000 failures occur
somewhere between 100 and 150 years. This accumulated
tectonic stress of
_
σ
0
T

t

t
0

1
:
0
1
:
5
MPa is small
compared with the simulated peaks in
σ
0
H

x

, which are of
the order 100
200 MPa.
_
σ
0
T

t

t
0

is also small compared
with the best-fitting
σ
0
B
that is presented in
The Size of
Stress at Different Scales
, of the order 60 MPa. Therefore,
_
σ
0
T

t

t
0

is typically less than 3% of
σ
0
B
for the duration
of the simulation.
In general, we assume that
σ
0
B
and
_
σ
0
T
may have differ-
ent orientations. For example, the principal compression of
the average background stress may be oriented nearly
perpendicular to the San Andreas fault (
Townend and
Zoback, 2004
); whereas the stress rate compression axis
can be approximately at a 45° angle (
Becker
et al.
, 2003
,
2005
) because shear on the San Andreas fault accommodates
most of the plate motion. Note, however, that the observation
of near perpendicular principal compression stress to the San
Andreas fault has been debated (
Hardebeck and Hauksson,
2001
,
2004
).
The term
σ
0
H

x

is the spatially varying deviatoric stress.
By definition, the expected value of each component in
σ
0
H

x

is zero. The heterogeneous stress is assumed to be due to all of
the stress changes caused by local inelastic deformations,
such as the slip distribution due to faulting, compaction,
fluids, thermal stresses, topography, and other factors. We as-
sume that the heterogeneity is described by two parameters:
1.
α
, where the amplitude spectrum of any
1D
cross section
through our
3D
σ
0
H

x

grid is proportional to
1
=k
α
r
, where
k
r

1
=r
is the wavenumber and
r
is a linear distance.
2. Heterogeneity ratio (
HR
) is a measure of the relative
amplitudes of the heterogeneous stress compared with
the uniform background stress. We measure the amplitude
of the stress tensor using the scalar inner product of the
deviatoric stress tensor with itself. This tensor inner pro-
duct (denoted by a colon) is equivalent to the second scalar
invariant of the deviatoric stress tensor,
I
0
2

σ
0
ij
:
σ
0
ij
(
Housner and Vreeland, 1965
).
I
0
2
is a measure of the shear
strain energy density.
HR


I
0
H
2
q

I
0
B
2
p



σ
0
H

x

:
σ
0
H

x

q

σ
0
B
:
σ
0
B
p
;
(2)
where

I
0
H
2
is the spatial average of the second invariant
of the heterogeneous deviatoric stress. Because
I
0
H
2
is
the sum of the squared components of
σ
0
H

x

, where
σ
0
H

x

σ
0

x
;
0

σ
0
B

σ
0

x
;
0


σ
0

x
;
0

, then

I
0
H
2
is the sum of the variances of the components of
σ
0

x
;
0

.
For this particular parameterization of the relative size of
the background stress (a constant) and the heterogeneous
stress (a spatially varying fractal-like term) to be independent
of the number of the points used in the simulation, the outer
scale must be fixed. For example, to achieve the same spec-
tral properties for a simulation of 1000 points length and a
simulation of 100,000 points length, we equate the length of
each simulation to the same outer scale. Then the spectral
properties up to a limiting spatial frequency will be the same.
This also means then that the distance measured between
grid points for the simulation with 100,000 points will be
100x smaller than the distance between grid points for the
simulation with 1000 points.
We assume the fractal-like statistics are stable for inter-
seismic times; therefore, the stochastic properties of
σ
0
H

x

,
described by
HR
and
α
, do not significantly evolve in time
for the simulations presented in the paper, and we do not
update
σ
0
H

x

after each event. Our focus is to compare
our results to stress inversions applied to background seismi-
city in between major seismic events over a time window in
the range of 1 to 150 years. The significant
3D
heterogeneous
stress changes that a major earthquake event introduces
would have to be taken into account to study the aftershock
period.
1400
D. E. Smith and T. H. Heaton
Assumptions/Limitations of This Stress Formulation
From the outset it is important to indicate clearly the
assumptions used in this paper and the possible limitations.
We do not attempt to create stress heterogeneity in
3D
from
first principals because of the inherent difficulties.
Aagaard
and Heaton (2008)
present numerical simulations of self-
sustaining heterogeneous slip and stress on a
2D
plane. Their
simulations were severely limited in spatial bandwidth by the
immense numerical calculations that are required to faithfully
create realistic
3D
stress heterogeneity. Furthermore, deriva-
tion of stress from dynamic modeling requires many assump-
tions, such as the distributions of fault orientations, fault
lengths, dynamic friction on faults, etc. Instead, we have
chosen to approach this problem with a simple statistical
model. On the positive side, this enables us to describe spa-
tially heterogeneous stress with two statistical parameters,
HR
and
α
, to generate synthetic focal mechanisms quickly and to
compare our simulations with real data to constrain the
statistical properties of the crust. On the other hand, this sta-
tistical approach makes many simplifying assumptions in an
attempt to obtain insight into the nature of stress variations in
the Earth
s crust and overlooks details that are necessary if one
wishes to model stress heterogeneity from first principles.
In particular, our spatial model of stress does not satisfy
the equations of static equilibrium. While we satisfy rota-
tional equilibrium (
σ
ij

σ
ji
), we do not satisfy the other
static equilibrium equation,
σ
ij;j

f
i
, where
f
i
is a body
force. In order to satisfy
σ
ij;j

f
i
, the solutions to which
are spatially smooth, and to have spatially heterogeneous
stress, we would have to include sources, which requires a
whole set of additional assumptions. Of course, it is the local
effect of sources that is the likely origin of heterogeneous
stress, but we are trying to avoid the complications of explic-
itly including these sources.
The other assumptions for the formulation presented in
this paper include: (1) there is no such thing as preexisting
planar faults in our model, which means that our seismicity
tends to cluster in
3D
clouds rather than along lineations or
planes, as seen in the real Earth. (2) When failure occurs at a
grid point, the stress only changes at that grid point (it
actually drops out of the simulation once it has failed). This
means that there is no explicit interaction between events.
This assumption is clearly inappropriate for large events that
change stress over a large region. However, our stress model,
given by equation
(1)
, is intended to simulate background
seismicity, where stress perturbations due to individual
events are small and should have little to no effect on the
other events included in the regional inversions. (3) We
assume failure occurs on fresh fracture, maximally oriented
planes at

45
° from the
σ
1
and
σ
3
principal stress axes. This
is a consequence of using a plastic yield criterion; however,
one would also expect fracture at

45
° from the principal
stresses in Coulomb stress if
μ

0
. In appendix c from
Smith (2006)
, a Coulomb failure criterion with nonzero fric-
tion and pressure is applied. Similar but more complicated
results are found compared with the plastic yield criterion
results. Specifically, failures are still biased toward the stres-
sing rate term, but there is more scatter because the conjugate
planes are no longer orientated 90° with respect to one
another. (4) Last, our stress heterogeneity is fractal-like; its
outer scale is set by the box size and the inner scale is set by
the resolution. Later, in the results found in
Comparison with
Hardebeck
s (2006) Analysis of Southern California Focal
Mechanisms
, we show that this outer scale in our simulations
can be matched to outer scales present in the Earth.
Creating the Fractal-Like, Heterogeneous
Stress,
σ
0
H

x

In creating the deviatoric heterogeneous stress term,
σ
0
H

x

, we start with a statistically isotropic stress, subtract
the small mean due to finite sample size, apply a spatial filter,
and then subtract the pressure. Namely, we use a Gaussian
random number generator to assign values at each grid point
to the six independent Cartesian components. This produces
Cartesian tensor components that can be described by a
Gaussian distribution and principal stresses that can be
described by a beta distribution. The standard deviations of
the diagonal terms in the tensor are assumed to be 1.0, and
the standard deviations of the off-diagonal terms are assumed
to be
1
=

2
p
; this scaling is required to produce isotropy (see
Data and Resources
for hypersphere point picking).
Now we specify the desired spatial correlation upon
application of the spatial filter. Because the heterogeneous
stress grid has no preferred coordinate frame, we can
uniquely specify its correlation properties by defining the
spectral properties along any line that bisects the grid. If
k
r

1
=r
is the wavenumber along any line in which
r


x
i
x
i
p
specifies distance along the line, then we want
the expected value of the amplitude spectrum along any line
to be
E
hj
^
σ
H
ij

k
r
ji  
1
j
k
r
j

α
i

j
E
hj
^
σ
H
ij

k
r
ji 
1

2
p

1
j
k
r
j

α
i
j:
(3)
Because each Cartesian tensor component is initially
generated according to a Gaussian distribution, it is simple
to create this desired spatial correlation along
1D
lines in
space. Take the Fourier transform of each tensor component,
multiply by a
1D
spatial wavenumber filter,
^
F

k
r

1
j
k
r
j

α
;
(4)
and then take the inverse Fourier transform of the result. This
filter leaves the zero-frequency properties unchanged, while
filtering the short wavelength variations according to a power
law. It creates a fractal-like distribution that is relatively free
of inherent length scales besides the outer scale of the
Models of Stochastic, Spatially Varying Stress in the Crust Compatible with Focal-Mechanism Data
1401
maximum grid size and inner scale of the grid resolution. See
Figure
4
for examples of filtered scalars along
1D
lines and
their associated amplitude spectrums.
In Figure
5
, the effect of filtering on
3D
stress tensor
orientations is shown along a
1D
line. The filtered stress ten-
sor orientations along a
1D
line are represented as rotations
relative to a reference stress tensor orientation. Because there
are two conjugate planes associated with each stress tensor
and two slip directions with each plane, there are a total of
four possible reference orientations and four possible final
orientations for each stress tensor. Given this, we plot all
16 possible rotations for each configuration pair. The
3D
rotations can be formulated several different ways. In this
figure we choose to represent each
3D
rotation as a single
rotation,
ω
, about a specified rotation axis,

θ
;
φ

.
θ
is the co-
latitude and
φ
is the longitude of the rotation axis in spherical
coordinates; therefore, in Figure
5
, the location of the circles
represents the rotation axis

θ
;
φ

, and the color of the circles
represents the amplitude of the rotation,
ω
. See the
Data and
Resources
section for a link to a description of how the an-
gles,
ω
,

θ
;
φ

, convert to strike, dip, and rake

Θ
;
δ
;
λ

.In
general,

ω
;

θ
;
φ

can be converted to a quaternion (a four
component vector that describes
3D
rotations) and then the
quaternion can be converted to

Θ
;
δ
;
λ

. The definition of
the components of the quaternion vector is
q
0

cos

ω
=
2

q
1

sin

ω
=
2

sin

θ

cos

φ

q
2

sin

ω
=
2

sin

θ

sin

φ

q
3

sin

ω
=
2

cos

θ

;
(5)
where
j
~
q
j

q
2
0

q
2
1

q
2
2

q
2
3
p

1
. The rotation in
terms of quaternions can then be translated into

Θ
;
δ
;
λ

as follows,
tan
Θ

q
0
q
1

q
2
q
3
q
0
q
2

q
1
q
3
tan
λ

q
0
q
1

q
2
q
3
q
0
q
2

q
1
q
3
tan
δ

2

q
0
q
1

q
2
q
3

=
sin
λ
q
0
q
0

q
1
q
1

q
2
q
2

q
3
q
3
:
(6)
To achieve the spectrum given by equation
(3)
for our
grid in
3D
, we use the following methodology and spatial
wavenumber filter. We take the Fourier transform of our
3D
spatial grid of random numbers,
σ
H
0
ij

x
1
;x
2
;x
3

, with
respect to the three Cartesian coordinates to obtain
^
σ
H
0
ij

k
1
;k
2
;k
3

. We then multiply by a
3D
spatial wavenum-
ber filter,
^
F

k
1
;k
2
;k
3

, and take the inverse Fourier trans-
form of the result to generate,
σ
H
ij

x
1
;x
2
;x
3

. For proper
scaling we also apply some normalization, where
x
i
is the
length scale considered,
x
0
i
is the minimum length scale
(the same for each dimension),
L
i
is the total length of
the grid in each dimension, and
N
i
is the number of points
in the grid in each dimension, so that
^
F

k
1
;k
2
;k
3


1


L
i
L
i
p

3
p

k
i
k
i
p


f

α
;n

3

where
k
i

1
x
i
and
L
i

N
i
x
0
i
:
(7)
The function
f

α
;n

3

is required to produce
1D
spectra with a power-law decay of
k

α
r
from a grid of
n
0
(a)
(b)
12345
x 10
6
–4
–2
0
2
4
Random Gaussian White Noise, Filtered with
α
= 0.0
Should Be Equivalent to No Filtering
Amplitude of a Scalar
Length
10
–7
10
–6
10
–5
10
–4
10
–3
Spatial Frequency
Fourier Amplitude of a Scalar
Random Gaussian White Noise, Filtered with
α
= 0.0
Should Be Equivalent to No Filtering
Random Gaussian White Noise, Filtered with
α
= 0.5
Random Gaussian White Noise, Filtered with
α
= 0.5
Amplitude of a Scalar
Length
Spatial Frequency
Fourier Amplitude of a Scalar
–4
–2
0
2
4
012345
x 10
6
10
–2
10
–1
10
–3
10
–2
10
–1
10
–3
10
–4
10
0
10
–7
10
–6
10
–5
10
–4
10
–3
Figure
4.
(a) Gaussian white noise in the top panel for a
1D
line with a spatial correlation parameter,
α

0
:
0
(equivalent to no filtering).
This could, for example, be one component of a heterogeneous stress tensor,
σ
H

x

. The bottom panel is a log-log plot of the Fourier
amplitude spectra of this noise versus its spatial frequency. (b) Gaussian white noise filtered with
α

0
:
5
in the top panel. The bottom
panel is the log-log plot of the Fourier amplitude spectrum of this filtered noise versus its spatial frequency. Note that the slope of the trend

0
:
5
. This is the desired
1D
slope for a spatial correlation parameter of
α

0
:
5
. The color version of this figure is available only in the
electronic edition.
1402
D. E. Smith and T. H. Heaton
dimensions (in our case,
n

3
). We find that the simplest
functional forms that approximate
f

α
;n

in both
2D
and
3D
are hyperbolic functions. By trial and error, aided by Kalei-
dagraph
s fitting function, we find the following hyperbolic
approximation for
f

α
;n

in
3D
:
f

α
;n

3

2
:
97


α

3
:
5
3
:
33

2

1

s

1
:
(8)
To employ the same type of operation in
2D
, we find a
2D
spatial wavenumber filter,
^
F

k
1
;k
2

, where
^
F

k
1
;k
2


1


L
i
L
i
p

2
p

k
i
k
i
p


f

α
;n

2

;
(9)
and
f

α
;n

2

1
:
53


α

1
:
87
1
:
57

2

1

s

1
:
(10)
Typically, these equations are only valid for
α
<
1
:
5
, because
at
α
1
:
5
the filtered Gaussian noise becomes so smooth
that the power-law linearity begins to break down on the
log spectral amplitude versus log spatial wavelength plots.
Figure
6
shows examples of
2D
cross sections through
3D
filtered stress. We apply equations
(4)
,
(7)
, and
(8)
to the
components of our stress tensors with initially random
Gaussian distributions in our
3D
grid. We then plot one com-
ponent,
σ
H
11

x

, along a
2D
cross section of our
3D
grid to
show how the filtering equations affect stress on a plane.
This filtered heterogeneous stress,
σ
H

x

, is constructed
so that the mean value of every component is zero, which
means that there are no preferred orientations of the stress.
The expected value of the unfiltered components (
α

0
)is
zero; however, for grids with a finite number of points, the
expected value of
j

R

x
j
goes as

2
p
=

π
N
p
, where
N
is the
number of grid points if
R

x

is a normal random variable
with a mean of
μ

0
:
0
and a standard deviation of
σ

1
:
0
.
When the filters in equations
(7)
or
(9)
are applied to the
unfiltered components, the amplitudes of the nonzero spatial
frequencies are reduced, resulting in a larger relative mean.
Measuring this nonzero mean as
Figure
5.
Plots of stress tensor orientations for
1D
lines, 100,001 points each with different spatial correlation parameters,
α
, applied.
Each stress tensor orientation is represented as a
3D
rotation relative to a reference orientation. The
3D
rotations are defined by three angles,
ω
,
which is the amplitude of the rotation and

θ
;
φ

which is the axis of rotation.

θ
;
φ

defines the spatial colatitude and longitude of the points in
the plots, and
ω
defines the color. Note that the color bar associated with
ω
ranges from 0 to
2
π
radians. Interestingly, for
α

0
:
0
, the
orientations and colors are fairly uniformly distributed. As the spatial correlation increases, the spatial clustering and color grouping of the
points increases, until for
α

1
:
5
, one can distinctly see lines on the unit sphere representing orientations. This is simply showing the fairly
smooth variation of stress tensor orientations along the modeled
1D
line when
α

1
:
5
. The color version of this figure is available only in the
electronic edition.
Models of Stochastic, Spatially Varying Stress in the Crust Compatible with Focal-Mechanism Data
1403
MeanBias



σ
0
H

x

:

σ
0
H

x

q
;
(11)
it is on the order of 1% or less for
0
<
α
0
:
8
and is on the
order of 5% for
α
1
:
0
. While these are small values, a
σ
H

x

with a zero mean is desired; therefore, we explicitly
subtract the small mean value of each component of
σ
H

x

prior to filtering so that the final, filtered
σ
H

x

has zero
mean components. Last, when
σ
H

x

has been filtered, we
subtract the pressure to create our deviatoric heterogeneous
stress
σ
0
H

x

for equation
(1)
.
Fracture Criterion Used to Create Synthetic
Focal-Mechanism Catalogs
To create synthetic failures and their associated earth-
quake focal mechanisms, we need to select an appropriate
fracture criterion. Our preferred fracture criterion is the
Hencky
Mises plastic yield condition (
Housner and Vree-
land, 1965
) because of its simplicity. It predicts failure when
I
0
2
, which is proportional to the shear strain energy, a scalar
quantity that is related to the maximum shear stress, is greater
than a threshold value. Because
I
0
2
and the maximum shear
stress are invariant quantities, this failure criterion works
regardless of the coordinate system or orientation of the
individual stress tensors. The coefficient of friction is essen-
tially zero (optimally oriented planes), and pressure does
not enter into the equation. (If one wishes to investigate non-
zero pressures and coefficients of friction see appendix c,
Coulomb Fracture Criterion, in
Smith, 2006
.) Last, because
we are dealing with optimally oriented planes, the conjugate
planes become mathematically indistinguishable.
According to
Housner and Vreeland (1965)
, failure
occurs when
I
0
2

x
;t

2
3
σ
2
0
;
(12)
where
σ
0
is the uniaxial yield stress and
I
0
2

x
;t

is the second
invariant of the deviatoric stress tensor,
σ
0

x
;t

, where
I
0
2

x
;t

σ
0

x
;t

:
σ
0

x
;t

:
(13)
Given the fracture criterion in equation
(12)
, definition of
I
0
2

x
;t

in equation
(13)
, and equation for deviatoric stress
tensor
(1)
, we show in Figures
7
and
8
the effect of our
two statistical parameters,
HR
and
α
, on generating failures
and their associated focal mechanisms.
Figure
6.
Plots of a scalar tensor component,
σ
H
11

x

, with different filters applied. The
2D
cross sections are
x
-
y
planes through ap-
proximately the center of the
3D
grid. We start with Gaussian noise and apply
α

0
:
0
, 0.7, 1.0, and 1.5 to produce scalar tensor components
with different spectral
1D
falloffs. The color version of this figure is available only in the electronic edition.
1404
D. E. Smith and T. H. Heaton