1

Thermodynamic Method for Generating Random 

Stress Distributions on an Earthquake Fault

By Michael Barall, Invisible Software, Inc., and Ruth A. Harris, U.S. Geological Survey

Open-File Report 2012–1226

U.S. Department of the Interior

U.S. Geological Survey
2

U.S. Department of the Interior

KEN SALAZAR, Secretary

U.S. Geological Survey

Marcia K. McNutt, Director

U.S. Geological Survey, Reston, Virginia: 2012

For more information on the USGS—the Federal source for science about the Earth,

its natural and living resources, natural hazards, and the environment—visit 

http://www.usgs.gov or call 1–888–ASK–USGS

For an overview of USGS information products, including maps, imagery, and publications,

visit http://www.usgs.gov/pubprod

To order this and other USGS information products, visit http://store.usgs.gov

Suggested citation:

Barall, Michael, and Harris, R.A., 2012, Thermodynamic method for generating random stress 

distributions on an earthquake fault: U.S. Geological Survey Open-File Report 2012–1226, 112 p.

Any use of trade, product, or firm names is for descriptive purposes only and does not imply 

endorsement by the U.S. Government.

Although this report is in the public domain, permission must be secured from the individual 

copyright owners to reproduce any copyrighted material contained within this report.
3

Contents

Abstract ......................................................................................................................................................... 6

Introduction.................................................................................................................................................... 7

Purpose and Scope ................................................................................................................................... 7

Previous Work............................................................................................................................................ 8

Thermodynamic Approach......................................................................................................................... 9

Thermodynamic Concepts........................................................................................................................... 11

Ising Model............................................................................................................................................... 11

Micro-States and Boltzmann Probability .................................................................................................. 11

Equipartition Theorem.............................................................................................................................. 12

Thermodynamic Analogy ............................................................................................................................. 14

Thermodynamic Fault Model ....................................................................................................................... 15

Basis Displacement Fields....................................................................................................................... 15

Properties of the Basis Displacement Fields............................................................................................ 16

Thermodynamic Model Construction........................................................................................................ 18

Ultraviolet Divergence.................................................................................................................................. 21

Blackbody Radiation ................................................................................................................................ 21

Finite Fault Strength Restricts the Set of Accessible States .................................................................... 22

Scaling Relations and Self-Similarity ........................................................................................................... 25

Scaling Relations ..................................................................................................................................... 25

Energy Rolloff and Power Spectrum........................................................................................................ 27

State Restriction and Probability Distribution............................................................................................... 28

Energy Rolloff Function............................................................................................................................ 28

State Restriction....................................................................................................................................... 28

Algorithm to Generate the Stress Field........................................................................................................ 30

Preliminary Steps..................................................................................................................................... 30

Generating Random Amplitudes and Phases .......................................................................................... 31

Summing the Random Stress Fields........................................................................................................ 32

Controlling Rupture Extent........................................................................................................................... 33

Depth Conditioning .................................................................................................................................. 33

Selection Process .................................................................................................................................... 34

Selection Algorithm.................................................................................................................................. 34

Correlation Length ................................................................................................................................... 36

Application to Smaller Earthquakes ......................................................................................................... 36

Post-Processing .......................................................................................................................................... 38

Nucleation.................................................................................................................................................... 39

The Problem of Nucleation....................................................................................................................... 39

The Day Radius ....................................................................................................................................... 40

Two-Stage Nucleation Method................................................................................................................. 40

Conclusions ................................................................................................................................................. 42

Acknowledgments ....................................................................................................................................... 43
4

Appendix A. Slider-Block Model .................................................................................................................. 47

Appendix B. Properties of the Basis Displacement Fields ........................................................................... 49

Fault Slip and Exponential Decay ............................................................................................................ 50

Derivatives of the Basis Displacement Fields .......................................................................................... 50

Elastostatic Equation ............................................................................................................................... 52

Elastic Energy .......................................................................................................................................... 52

Orthogonality............................................................................................................................................ 54

Normal Stress .......................................................................................................................................... 55

Shear Stress ............................................................................................................................................ 56

Total Shear Stress and Completeness..................................................................................................... 58

Appendix C. Derivation of the Earthquake Rate Formula ............................................................................ 60

Appendix D. Derivation of the Power Spectral Density ................................................................................ 62

Notations and Definitions ......................................................................................................................... 62

Filtered Stress Fields ............................................................................................................................... 63

Expected One-Point Statistics of Filtered Stress...................................................................................... 63

Small-Scale Stress Perturbations ............................................................................................................ 64

Method of Calculation .............................................................................................................................. 65

Calculation Assuming Gaussian Distribution............................................................................................ 66

Calculation Assuming Laplace Distribution .............................................................................................. 68

Calculation Assuming Exponential Distribution ........................................................................................ 69

Calculation Assuming Pareto Distribution ................................................................................................ 70

Technical Issues for One-Point Statistics................................................................................................. 71

Removing Infinities ............................................................................................................................... 71

Stress Distributions and Wave Amplitudes........................................................................................... 72

Fault Surface and Mean Values ........................................................................................................... 73

Product Space...................................................................................................................................... 74

Definition of Expected One-Point Statistics .......................................................................................... 74

Mean and Variance of Expected One-Point Statistics .......................................................................... 76

Appendix E. Algorithms for Energy Rolloff................................................................................................... 78

Energy Rolloff Function............................................................................................................................ 78

State Restriction....................................................................................................................................... 79

Appendix F. Mathematics of the Selection Algorithm................................................................................... 83

The Convolution....................................................................................................................................... 83

The Phase Shift ....................................................................................................................................... 84

Appendix G. Post-Processing...................................................................................................................... 85

Stress Field Rotation................................................................................................................................ 85

Filtering .................................................................................................................................................... 85

Scaling ..................................................................................................................................................... 86

Hypocenter Selection............................................................................................................................... 87

Size of Nucleation Zone........................................................................................................................... 87

Depth Conditioning .................................................................................................................................. 88

Scaling by Normal Stress......................................................................................................................... 89

Example of Scaling .................................................................................................................................. 89

Appendix H. One-Point Statistics................................................................................................................. 92

Discussion................................................................................................................................................ 92
5

Figures

Figure 1. Components of dynamic and kinematic earthquake simulations. .............................................. 100

Figure 2. Two-dimensional Ising model .................................................................................................... 101

Figure 3. Three-dimensional coordinate system....................................................................................... 102

Figure 4. Blackbody radiation. .................................................................................................................. 103

Figure 5. Energy rolloff function................................................................................................................ 105

Figure 6. Examples of random stress distributions generated by our method. ......................................... 107

Figure 7. Selection process for controlling rupture extent along-strike. .................................................... 108

Figure 8. Selection template. .................................................................................................................... 109

Figure 9. Effect of varying the correlation length....................................................................................... 111

Figure A–1. Slider-block model of an earthquake fault. ............................................................................ 112

Tables

Table 1. Mathematical notations............................................................................................................ 93

Table 2. Summary of thermodynamic analogies.................................................................................... 97

Conversion Factors

SI to Inch/Pound

Multiply By To obtain

Length

meter (m) 3.281 foot (ft) 

kilometer (km) 0.6214 mile (mi)

Pressure

megapascal (MPa) 145.0 pound per square inch (lb/in2)
6

Thermodynamic Method for Generating Random 

Stress Distributions on an Earthquake Fault

Michael Barall, Invisible Software, Inc., and Ruth A. Harris, U.S. Geological Survey

Abstract

This report presents a new method for generating random stress distributions on an 

earthquake fault, suitable for use as initial conditions in a dynamic rupture simulation. The 

method employs concepts from thermodynamics and statistical mechanics. A pattern of fault slip 

is considered to be analogous to a micro-state of a thermodynamic system. The energy of the 

micro-state is taken to be the elastic energy stored in the surrounding medium. Then, the 

Boltzmann distribution gives the probability of a given pattern of fault slip and stress. We show 

how to decompose the system into independent degrees of freedom, which makes it 

computationally feasible to select a random state. However, due to the equipartition theorem, 

straightforward application of the Boltzmann distribution leads to a divergence which predicts 

infinite stress. To avoid equipartition, we show that the finite strength of the fault acts to restrict 

the possible states of the system. By analyzing a set of earthquake scaling relations, we derive a 

new formula for the expected power spectral density of the stress distribution, which allows us to 

construct a computer algorithm free of infinities. We then present a new technique for 

controlling the extent of the rupture by generating a random stress distribution thousands of 

times larger than the fault surface, and selecting a portion which, by chance, has a positive stress 

perturbation of the desired size. Finally, we present a new two-stage nucleation method that 

combines a small zone of forced rupture with a larger zone of reduced fracture energy.
7

Introduction

Dynamic earthquake simulations are powerful computational tools. They incorporate our 

knowledge of earthquake physics and allow for the prediction of earthquake magnitudes and 

effects. These simulations also require many assumptions about how faults and earthquakes 

work, including values for the initial stress conditions, and parameters that describe the fault 

geometry, the frictional failure criterion, and the material properties (for example, Harris (2004)).

Fault geometry and material properties can, in principle, be estimated given adequate 

geophysical measurements, but other aspects of dynamic earthquake rupture remain more 

elusive. The most challenging parameters are those that relate to fault friction and the initial 

stress conditions on faults. We leave fault friction for others to study and refer the reader to 

overviews by Tullis (2007), Lapusta (2009), and Bizzarri (2011).

In this report we concentrate on the problem of defining appropriate probability 

distributions for fault stress heterogeneity. We present a method for prescribing random 

realizations of stress at the onset of dynamic earthquake rupture simulations. Preliminary tests of 

this approach show the resulting heterogeneous initial-stress distributions to operate well in 

computer codes that are included in the Southern California Earthquake Center (SCEC) Dynamic 

Rupture Code group (http://scecdata.usc.edu/cvws; Harris, Barall, Archuleta, and others (2009);

Harris, Barall, Andrews, and others (2011)). Computer code to implement the method 

accompanies this report.

Purpose and Scope

We are concerned with setting up the initial conditions for a dynamic rupture simulation. 

Detailed examples of initial conditions for dynamic rupture simulations are available on the 

SCEC website (http://scecdata.usc.edu/cvws). To perform a dynamic rupture simulation, one 

must specify:

 Fault geometry.

 Rock properties (which we presume to be linear elastic).

 Frictional constitutive properties of the fault (which may have parameters that vary 

with position on the fault surface).

 Initial normal stress at each point on the fault surface.

 Initial shear stress at each point on the fault surface.

In this work, we give a method for randomly generating the initial shear stresses. All the 

other components listed above are presumed to be specified in advance.

Our method addresses the case of a planar fault, embedded in a uniform linear elastic 

half-space, with uniform orientations of the initial shear stress. For ease of discussion, we 

describe the method for the case of a vertical planar fault with horizontal initial shear stresses. 

There is nothing in the development of the method that relies on the fault being vertical or the 

shear stress orientation being horizontal, and so the method can be applied unchanged to dipping 

faults and to faults with non-horizontal (uniform) initial shear stress orientations.

For ease of discussion, we also assume that the simulated earthquake is large enough to 

rupture the full depth of the seismogenic zone and that the fault obeys a linear slip-weakening 

friction law. For the case of earthquakes too small to rupture the full depth of the seismogenic 
8

zone, a minor modification to the method is required, which we describe. The method does not 

depend on details of the friction law and so may be applied unchanged to other friction laws.

We do not address the cases of non-planar faults, non-uniform initial shear stress 

orientations, or non-uniform linear elastic rock properties. One can of course apply our method 

(or any other method) to such cases simply by generating an initial shear stress distribution for a 

similar case where the method is applicable, and then just using that stress distribution despite 

the fact that the method employed is not directly applicable to the case at hand. In doing so, 

however, one must be aware that the method is being employed in a situation where it is not 

specifically applicable.

Our method cannot be used for plastic rock properties, because that would require a 

randomly generated initial stress tensor throughout the three-dimensional medium. Our method 

only generates a random initial shear stress distribution on the fault surface.

Previous Work

Heterogeneous fault-stress distributions and heterogeneous fault-slip distributions have 

been studied for at least 35 years, including early work by Yamashita (1976), Das and Aki 

(1977), Mikumo and Miyatake (1978), Hanks (1979), McGarr and others (1979), Andrews 

(1980; 1981), Boatwright (1981), and Day (1982). More recent studies include those by Beroza 

and Mikumo (1996), Bouchon (1997), Day and others (1998), Mai and Beroza (2002), Oglesby 

and Day (2002), Guatteri and others (2003), Guatteri, Mai, and Beroza (2004), Lavallée and 

others (2006), Ripperger and others (2007), Ripperger, Mai, and Ampuero (2008), Aagaard and 

Heaton (2008), Schmedes et al. (2010), and Andrews and Barall (2011), among others. Some of 

these studies have parameterized the fault slip (or slip velocity), others the fault stress. The 

former is needed for kinematic models of earthquake rupture, in which the timing and amount of 

slip at each point on a fault are pre-determined before the start of the simulation. The latter is 

required for spontaneous models of dynamic earthquake rupture, in which the rupture extent and 

the timing and amount of slip are determined by the physics of the processes occurring during 

the course of the simulation (fig. 1). Among the authors who have worked with heterogeneous 

initial stresses in spontaneous rupture simulations, the need for heterogeneity has been explained 

as necessary to produce reasonable rupture propagation direction and fault slip patterns (for 

example, Yamashita, 1976), and as an explanation for variable rupture speeds and spatial 

variations of peak slip velocity (for example, Day, 1982). Some authors have proposed that the 

heterogeneity may reside in the initial, static, or dynamic frictional stresses or all three ( for 

example, Mikumo and Miyatake, 1978; Peyrat and others, 2001; Bouchon, 1997; Mavroeidis and 

Papageorgiou, 2010), with Mikumo and Miyatake (1978) proposing that the heterogeneity may 

be due to geometrical fault complexity and may explain observations of high-frequency radiation 

patterns and locations of aftershocks. Many have mentioned that heterogeneity in stresses is

essential for multi-earthquake cycle simulations.

This report does not assume that we know the exact cause of spatial fluctuations in the 

initial stresses on faults, it just recognizes that realistic ground motions are unlikely to be 

produced by assuming constant slip or constant initial stress conditions on a planar fault surface.

We therefore provide a formulation to describe the heterogeneity. As we improve our 

understanding of the fine-scale details of faults, along with our ability to realistically simulate 

them ( for example, Kaneko and others, 2010; Shi and Day, 2011), the approach presented here 

may evolve. For now, however, we need to adopt a stochastic method to generate heterogeneities 

in the initial stress state in dynamic, spontaneous earthquake rupture simulations.
9

Thermodynamic Approach

Many existing methods for generating random stress distributions attempt to emulate 

statistical properties of earthquakes inferred from slip inversion models (for example, Mai and 

Beroza, 2002; Lavallée and others 2006), or from considerations of earthquake self-similarity 

(for example, Andrews and Barall, 2011). Here we explore a new approach, which is to generate 

random stress distributions using concepts and tools drawn from thermodynamics and statistical 

mechanics. Our approach is motivated by the two-dimensional Ising model, which is an example 

of how a physical system with simple laws can self-organize into large structures. Smalley and 

others (1985) and Huang and others (1998) construct models of stress transfer on an earthquake 

fault with features analogous to the Ising model. Other authors who have applied statistical 

mechanics concepts to the study of earthquake faults include Rundle and others (2003) and 

Abaimov (2009). Here we employ an analogy to the Ising model to construct a method for 

generating random heterogeneous initial stress.

In statistical mechanics, it is presumed that the precise state of a physical system is 

unknowable, so one must talk about a probability distribution of states. The probability that a 

system is in a given state depends on the state’s energy; the higher the energy, the lower the 

probability. We consider a pattern of slip on a fault surface to be analogous to a thermodynamic 

state. Fault slip distorts the three-dimensional elastic medium surrounding the fault, and we 

consider the resulting elastic energy to be analogous to the energy of a thermodynamic state. The 

goal of this report is to carry this analogy through to obtain a probability distribution for the 

resulting patterns of stress on the fault.

Our work begins with a review of the two-dimensional Ising model and the 

thermodynamic concepts we need. Then, we introduce a set of three-dimensional displacement 

fields, which provide a mathematically tractable way to describe patterns of fault slip and stress, 

and the resulting distortion of the three-dimensional elastic medium. The displacement fields 

have the property of orthogonality, which means that when combined in a linear superposition 

their elastic energies are additive. So, these displacement fields are analogous to the degrees of 

freedom of a thermodynamic system, and we can use the equipartition theorem to develop a 

probability distribution.

Unfortunately, the equipartition theorem predicts infinitely large stresses due to the everincreasing contributions of short-wavelength displacement fields. This phenomenon, called 

ultraviolet divergence, occurs in some statistical mechanics problems such as the calculation of 

the spectrum of blackbody radiation. Its resolution is to restrict the possible states of the system. 

For us, the finite strength of the fault restricts the possible patterns of fault slip and stress, 

because a stress pattern that exceeds the yield strength of the fault cannot occur. By careful 

consideration of earthquake scaling relations and self-similarity, we derive a formula showing 

how the finite fault strength causes the energy distribution to deviate from equipartition. This 

allows us to develop a probability distribution that predicts finite levels of stress. We then 

present an algorithm for generating random realizations of initial stress according to this 

probability distribution.

After generating a random stress field, several post-processing steps are required to make 

it suitable for use in a dynamic rupture simulation. It is desirable to make the simulated rupture 

stop spontaneously by encountering low stress outside the intended rupture area. We achieve this 

using a new method, which generates a random stress field over an extremely large area, and 

then selects a small portion of the random stress field to use on the simulated fault surface,

which, just by chance, has the required low-stress regions. Then, we indicate how to scale and 
10

filter the stress field, to make it consistent with the assumed fault friction law, and to suppress 

stress fluctuations too small to be resolved by the dynamic rupture code. Finally, we describe a 

new technique for initiating, or nucleating, the simulated rupture, which combines a small region 

of forced rupture with a larger region of reduced fracture energy.
11

Thermodynamic Concepts

Our method is based on concepts from thermodynamics, and in particular the twodimensional Ising model. In this section we review the relevant concepts. For a general reference 

see Reif (1965).

There are several reasons why we have chosen to use thermodynamics as a basis for our 

method. First, the Ising model is a classic example of how a stochastic physical system can selforganize into large structures. For the Ising model, the large structures are magnetic domains in a 

ferromagnetic material. For an earthquake fault, the large structures are stress perturbations large 

enough to sustain a big earthquake. Second, the principles of thermodynamics are broadly 

applicable to many physical systems, and so we may expect this approach to yield physically 

reasonable results. And third, this approach lets us utilize known results and concepts from 

thermodynamic, such as the Boltzmann probability distribution and the equipartition theorem, in 

the construction of our method.

Ising Model

The two-dimensional Ising model (Ising, 1925; Onsager, 1944) is a simple model of a 

magnetic material, such as iron (fig. 2). Each point on a two-dimensional lattice holds a spin that 

can point in one of two directions, either up or down. Each spin represents, for example, a single 

atom of iron which behaves like a tiny magnet. The north pole of the magnet can point either up 

or down. (Quantum mechanics limits the orientation of the spin to one of two discrete 

directions.) The value of the nth spin is denoted by , which is if the spin points up or if 

the spin points down.

In the Ising model, each spin interacts only with its four nearest neighbors. (If the lattice 

is of finite size, then spins on the edge of the lattice have only three or two nearest neighbors. 

Alternatively, one can impose periodic boundary conditions, in which case every spin on a finite 

lattice can have four nearest neighbors.) The total energy, , of the system is then given by the 

formula

 ∑ 

 

(1)

where the sum runs over nearest-neighbor pairs. For a ferromagnetic material, the coupling 

constant, is positive, which makes it energetically favorable for adjacent spins to line up. (For 

a paramagnetic material, is negative.) It can be shown that if the temperature is not too high, 

then the spins in the Ising model spontaneously organize themselves into large domains in which 

most of the spins are aligned, which in turn can create a macroscopic magnetic field.

Micro-States and Boltzmann Probability

A micro-state of a system is a complete specification of the system’s microscopic state. 

For example, a micro-state of the Ising model would specify the values of all the spins . In 

general, one cannot know the micro-state of a system. Instead, one knows macroscopic 

parameters such as the temperature. Given the values of macroscopic parameters, one then infers 

a probability distribution for the possible micro-states.
12

The Boltzmann probability distribution gives the probability that the system is in 

micro-state :

 

 

 

 

 here 

 

 

 (2)

Here is the energy of the system when it is in micro-state , and is proportional to the 

reciprocal of the temperature, . The normalization factor is called the partition function

and is defined to be

 ∑ 

 

 

(3)

where the sum runs over all accessible micro-states . If the micro-states form a continuum, then 

the sum is replaced with an integral. A useful feature of the partition function is that the expected 

energy of the system, denoted 〈 〉, can be computed by differentiating the partition function:

〈 〉 

 

 ∑ 

 

 

 

 

 

 

 (4)

The Boltzmann distribution can be interpreted as expressing two principles: (1) physical 

systems tend to assume lower-energy states; and (2) if you pump a certain amount of energy into 

a physical system (for example, by heating it), then all accessible micro-states with that amount 

of energy, or less, are about equally likely.

Equipartition Theorem

The equipartition theorem is a consequence of the Boltzmann distribution (Boltzmann, 

1876). Suppose that (1) the system can be decomposed into a set of degrees of freedom, and (2) 

the total energy of the system is equal to the sum of the energies associated with each degree of 

freedom. Then, the Boltzmann probability distribution applies separately to each degree of 

freedom, as if they were independent random variables. So, if is a decomposition of the 

micro-state into a set of degrees of freedom, and if the energy of the micro-state can be expressed 

in the form

 ∑ 

 

(5)

then the probability of finding the system in micro-state is

 ∏ 

 

 here 

 

 

 

 

 

 

 (6)

This equation is the formula for the joint probability distribution of independent random 

variables.
13

Continuing the equipartition theorem, in addition suppose that (3) each degree of freedom 

is a continuous variable that can assume any real value, and (4) the energy associated with an 

individual degree of freedom depends quadratically on its value. Then, each degree of freedom 

has the same expected energy,

〈 〉 

 

 

for all (7)

where is the energy associated with the nth degree of freedom.

Other thermodynamic concepts, including density of states and ultraviolet divergence, 

will be introduced as needed.
14

Thermodynamic Analogy

Here we give an overview of how we intend to construct a heterogeneous stress pattern 

by analogy to the Ising model. In the next section we begin to develop the machinery to carry out 

this plan. The steps to construct the thermodynamic analogy are:

1. The system is the earthquake fault and surrounding rock.

2. A micro-state is defined to be a pattern of slip on the fault.

3. The energy, , of a given micro-state is defined to be the elastic energy associated 

with the given pattern of slip.

4. Choose a thermodynamic parameter, (This will turn out to be equivalent to 

choosing a correlation length for the pattern).

5. Use a computer to randomly construct a micro-state, so that the probability of 

choosing a given pattern of slip equals its Boltzmann probability (which depends on 

its energy).

6. To make the computation feasible, decompose the system into a set of degrees of 

freedom, so that the amplitude for each degree of freedom can be selected 

independently according to the Boltzmann distribution.

7. Compute the pattern of shear stress resulting from the selected pattern of slip.

Table 2 summarizes the analogy between a thermodynamic system and an earthquake 

fault.

In attempting to carry out this plan, we encounter a major obstacle: the Boltzmann 

distribution predicts infinitely large stresses. This behavior is known as an ultraviolet divergence.

To make the stresses finite, we need to revisit the Gutenberg-Richter relation and the subject of 

earthquake self-similarity (or scale-invariance), and derive a new result concerning the stress 

distribution that is implied by self-similarity.

Once the ultraviolet divergence is resolved, the above steps can be used to generate a 

stress distribution. Then, we show how to control the extent of the rupture, select the correlation 

length, and nucleate the rupture.
15

Thermodynamic Fault Model

We must work in three dimensions to construct a thermodynamic fault model with 

correct energy spectrum. (Appendix A explains why a two-dimensional fault model cannot be 

used.) To work in three dimensions, we introduce a set of three-dimensional displacement fields, 

which are approximate solutions to the problem of fault slip in an elastic medium. We can then 

describe a three-dimensional system configuration by forming a linear combination of the 

displacement fields. These fields are orthogonal, which means that the elastic energy associated 

with a linear combination of displacement fields equals the sum of the elastic energies of the 

individual fields. Due to their orthogonality, these fields act as the degrees of freedom of a 

thermodynamic system.

Basis Displacement Fields

The setting is a vertical, strike-slip fault in a three-dimensional elastic full-space (fig. 3). 

Points in space are denoted by coordinates where the coordinate is vertical, the 

coordinate is parallel to the fault strike, and the coordinate is perpendicular to the fault. The 

fault plane is .

We impose periodic boundary conditions in the and directions (but not in the 

direction). Let and be the periods in the and directions, respectively. Then the fault 

area is .

A displacement is represented by a vector . Our basis displacement fields are 

defined as:

 

   

 

 

 

 

 

 | |

 ( 

 )  (8)

 

   

 | |

 ( 

 )  (9)

 

   

 

 

 

 

 

 | |

 ( 

 ) (10)

Here 

 is an arbitrary amplitude,  is an arbitrary phase, and are 

wavenumbers in the and directions, and is the scalar wavenumber defined as:

 √ 

  

 (11)

To understand the structure of these basis displacement fields, take a closer look at the 

formula for , which is displacement in the along-strike direction. It is a product of four terms:

 An arbitrary amplitude 

 .

 A term 

 | | which decays exponentially with increasing distance from the fault, 

over a characteristic distance .

 A term ( 

 ) which is a sinusoidal variation in the fault 

plane, with wavelength .
16

 A term 

 which changes sign when crossing the fault plane , producing 

a discontinuity at the fault plane which represents fault slip.

The formulas for and are needed to satisfy the elastostatic equation and to make the 

shear stress point along-strike. This is discussed further below.

Properties of the Basis Displacement Fields

We now discuss the properties of the basis displacement fields. Mathematical proofs are 

in appendix B. Mathematical notations are summarized in table 1.

Property 1. Away from the fault, the displacements are an exact solution of the 

elastostatic equation:

 (12)

 (13)

 

 

 

( ) (14)

Here is the stress tensor, is the strain tensor, is the shear modulus, and is the 

volumetric Lame modulus.

Property 2. The shear stress, , is the same on both sides of the fault, as it should be.

Property 3. The normal stress, , is not the same on both sides of the fault. This is what 

makes the displacement field an approximate, rather than exact, solution to the elastostatic fault 

slip problem.

Property 4. At the fault, the displacements describe fault slip which varies sinusoidally, 

with wavenumbers and in the along-dip and along-strike directions, and with amplitude 

proportional to .

Property 5. Away from the fault, the displacements die off exponentially over 

characteristic distance .

This shows that fault slip with long wavelength (small ) penetrates more deeply into the 

three-dimensional medium than fault slip with short wavelength (large ).

Property 6. The elastic energy in the three-dimensional medium is:

 

 

 

 

 (

 

 

 

 

 

)

 

 

 

 

 (15)

This is the expected energy spectrum, with 

 

.

Property 7. The displacement fields are orthogonal. If you form a linear combination of 

these displacement fields, with different wave vectors, then the total elastic energy, , of the 

linear combination equals the sum of the energies of the individual fields:

 total 

 

 

 ∫ (16)
17

Orthogonality makes the basis displacement fields satisfy the premises of the 

equipartition theorem. This is discussed below.

Remark: At this point we pause to explain the meaning of an integral over wavenumber 

space, such as the one shown above. Because we impose periodic boundary conditions, the 

wavenumbers are restricted to a set of discrete values

 

 

 

and 

 

 

(17)

where and are integers. An integral over wavenumbers and is defined to be a sum 

over the admissible discrete values:

∫ 

 

 

 

 

∫ 

 

 

 

 

 

∑ 

 

 

(18)

The factor 

 

is called the density of states because it is the number of admissible 

wavenumber values per unit area in the -plane. Density of states is a standard concept from 

thermodynamics.

The wave vectors and describe the same physical state. To see this, 

note that the values of , , and are unchanged if we make the replacements and 

 , provided that we take 

  and   . 

Therefore, when integrating over wavenumber space, the domain of integration includes only 

one-half of the -plane, chosen so that only one member of each pair, and 

 , lies in the domain of integration. Also, the domain of integration excludes the zero 

wave vector 

 , because the zero wave vector contributes zero energy and zero 

stress.

If the value of a function does not depend on the direction of the wave vector, 

then is defined to be isotropic. In this case, we can write 

 , and we have

∫ 

 ∫ 

 

 

 (19)

The coefficient on the right-hand side is , rather than , because we only integrate over onehalf of the -plane.

Property 8. The shear stress on the fault is parallel to the direction, that is, it is oriented 

along-strike. The shear stress is:

 

 √

 

 

 

 ( 

 )  (20)
18

The above formula shows that the set of basis displacement fields is complete in the 

following sense: Given any pattern of shear stress oriented along-strike, there is a linear 

combination of basis displacement fields which produces that pattern of shear stress. (This 

statement assumes that the pattern of shear stress is square-integrable, that is, the pattern of 

shear stress has a finite mean-square value.) In fact, the linear combination of basis displacement 

fields is effectively a Fourier transform of the shear stress pattern.

Property 9. The mean-square shear stress is:

 rms

 

  

 

 

 

 (21)

This equation is just the mean-square of the previous equation. The mean-square shear stress, 

 

 

  , is called the power spectral density of the stress pattern.

Property 10. If you form a linear combination of these displacement fields, with different 

wave vectors, then the mean-square shear stress, 

 

, of the linear combination equals the sum 

of the mean-square shear stresses of the individual fields:

 total

 

 

 

 ∫ rms

 

 

 

 

 ∫ 

 (22)

This property is just a statement of Parseval’s theorem.

Thermodynamic Model Construction

We now show how to use the basis displacement fields to construct a thermodynamic 

fault model.

A micro-state of the thermodynamic system is a pattern of slip on the fault surface, and 

the corresponding displacement of the surrounding elastic medium, under the constraint that the 

shear stress must be oriented along-strike. We describe a micro-state as a linear combination of 

the basis displacement fields. In other words, we specify a micro-state of the system by giving 

the amplitude, 

 , and the phase,  , for each wave vector, .

We noted earlier that the wave vectors and describe the same 

physical state. So, when we form a linear combination of the basis displacement fields, we only 

include one member from each such pair of wave vectors. We also exclude the zero wave vector 

  because it describes a state with zero energy and zero stress.

Because each displacement field contains both an amplitude, , and a phase, , it 

represents two degrees of freedom in the thermodynamic system. This can be seen by 

considering the trigonometric identity

 

  (23)

where and . This shows that each basis displacement field can be 

expressed as a linear combination of two displacement fields with arbitrary amplitudes and 

 , which means that there are two degrees of freedom.

We now show that the equipartition theorem applies to our system. Recall the four 

premises of the equipartition theorem:
19

1. The system can be decomposed into a set of degrees of freedom.

2. The total energy of the system is equal to the sum of the energies associated with each 

degree of freedom.

3. Each degree of freedom is a continuous real variable which can assume an arbitrary 

value.

4. The energy of each degree of freedom depends quadratically on its value.

Premise 1 is satisfied because a micro-state is described by a sum of basis displacement 

fields. The completeness property of the basis displacement fields (property 8) states that any 

possible pattern of fault slip oriented along-strike can be produced by some linear combination of 

the basis displacement fields. This shows that our basis displacement fields provide a complete 

description that includes every possible micro-state.

Premise 2 is satisfied because of the orthogonality property of the basis displacement 

fields (property 7).

Premise 3 is satisfied because the amplitude, , of a basis displacement field can assume 

any real value, and the phase can assume any value from to . (This is equivalent to saying 

that the amplitudes and mentioned in the equation above can each assume any real value.)

Premise 4 is satisfied because the energy of a basis displacement field (property 6) is 

proportional to 

 

. (This is equivalent to saying that the two displacement fields on the righthand side of the equation above contribute energy proportional to 

 

and 

 

respectively, 

because 

 

  

.)

The first conclusion of the equipartition theorem says that the Boltzmann distribution 

applies separately to each degree of freedom, as if the degrees of freedom are independent 

random variables. We use this in the following way. Our goal is to randomly generate a pattern 

of shear stress on the fault surface, by selecting a micro-state chosen at random according to the 

Boltzmann probability distribution. The equipartition theorem states that we can make the 

random selection by independently choosing an amplitude, 

 and phase, , for 

each wave vector, , according to the Boltzmann distribution. The micro-state is the sum 

of the resulting basis displacement fields, and we can then compute the pattern of shear stress 

using the stress formula (property 8) and energy formula (property 6).

The ability to select amplitudes and phases independently for each wave vector is critical, 

because it makes our method computationally feasible.

(Remark: In the computer implementation of the method, we randomly choose an energy,

 , and phase, , for each wave vector; it is not necessary to compute the 

amplitude, .)

The second conclusion of the equipartition theorem says that each degree of freedom has 

the same expected energy of , where is the thermodynamic parameter that plays the role 

of reciprocal temperature. Because each of our basis displacement fields represents two degrees 

of freedom, each basis displacement field has expected energy equal to . That is,

⟨  ⟩ 

 

 

for all 

(24)

where the angle brackets indicated expected value.

This is a problem. To see why, use property 10 to write the expected value of the meansquare shear stress:
20

⟨ total

 

⟩ 

 

 

 ∫⟨ 

 ⟩ (25)

If 〈 〉 is a constant, then the integral is strongly divergent, which means that the meansquare shear stress is infinite. So, the Boltzmann distribution predicts an infinite amount of 

stress, which is clearly not acceptable.

This is called an ultraviolet divergence because the divergence occurs in the limit of large 

 or short wavelength. To proceed, we must find a way to resolve the ultraviolet divergence.
21

Ultraviolet Divergence

The ultraviolet divergence that we encountered is not a new problem. In some prior 

methods for generating heterogeneous initial stress conditions, the power spectral density 

behaves asymptotically like a power of the wavenumber in the limit of large :

 rms

 

  

 

 

 

as (26)

Any such method has an ultraviolet divergence if , because in that case 

 

falls 

off so slowly that it is not integrable over the -plane. For example, Andrews and Barall 

(2011) describe a method with . They resolve the divergence simply by imposing a highfrequency cutoff, limiting to a maximum value chosen so that the shortest wavelengths 

considered are several times the resolution of the dynamic rupture simulation code.

We consider an arbitrary high-frequency cutoff to be incompatible with our 

thermodynamic method. Such a cutoff is not motivated by any thermodynamic principle, results 

in a stress state that is not well-defined independently of the discretization, and as a practical 

matter places too much of the stress variability into the shortest allowed wavelengths (that is, a 

large fraction of the power in the stress variation is concentrated in wavelengths close to the 

smallest wavelengths that can be resolved by the simulation code). So, we seek to resolve the 

ultraviolet divergence without the imposition of a cutoff.

In this section, we present a well-known example of a thermodynamic system that 

exhibits an ultraviolet divergence—blackbody radiation. We show that the divergence can be 

resolved by restricting the set of micro-states of the system. Then, we show in a general way how 

the same concept can be applied to an earthquake fault, with the finite strength of the fault acting 

to restrict the set of micro-states.

Blackbody Radiation

The term ultraviolet divergence (sometimes called ultraviolet catastrophe) gets its name 

from a problem in statistical mechanics: computing the spectrum of blackbody radiation.

The setting is a container which is maintained at a temperature, (fig. 4). A small 

viewport allows the experimenter to observe the electromagnetic radiation inside the container. 

The problem is to calculate the amount of electromagnetic energy at each wavelength, and the 

total amount of electromagnetic energy.

A micro-state of the electromagnetic field can be described as a linear superposition of 

plane waves. A plane wave is characterized by a wave vector, polarization, and amplitude. The 

total energy of the electromagnetic field equals the sum of the energies of the individual plane 

waves. It can be shown that the electromagnetic field has two degrees of freedom for each wave 

vector and polarization (for example, Reif, 1965). So according to the equipartition theorem, the 

expected energy in the plane wave with wave vector and polarization is:

⟨ ⟩ 

 

 

 (27)

This immediately implies an ultraviolet divergence. Because there are an infinite number 

of admissible wave vectors, the equipartition theorem predicts that the total electromagnetic 
22

energy is infinite. Also, because the number of admissible wave vectors increases as the 

wavelength decreases, it predicts that the energy increases without limit as the wavelength 

decreases.

The ultraviolet divergence is resolved by introducing quantum mechanics, as shown by 

Planck. The key idea is to reduce the set of accessible micro-states. Instead of allowing the 

energy to have an arbitrary real value, quantum mechanics restricts the energy to a set of discrete 

values, called quantum levels:

 (28)

where is an integer, is Planck’s constant, is the speed of light, and is the scalar 

wavenumber. (Another way of saying this is that each plane wave contains an integral number of 

photons, where each photon carries energy, .)

It is important to note that the probability of a micro-state is still given by the Boltzmann 

probability distribution, even when quantum mechanics is introduced. The difference is that now 

the Boltzmann probability distribution is applied over a much smaller set of micro-states. With 

this change, the expected energy is given by the Planck distribution:

⟨ ⟩ 

 

 

 

 

 

 (29)

This resolves the ultraviolet divergence because the expression on the right-hand side has 

a finite sum, so the total energy is finite.

Let us look more closely at the expression 

 appearing on the righthand side. For small or long wavelength, it approaches , which means that the expected 

energy is close to classical value . But for large or short wavelength, it approaches , 

which means that the expected energy declines to zero.

Finite Fault Strength Restricts the Set of Accessible States

Our objective is to resolve the ultraviolet divergence for earthquake faults in a manner 

analogous to the way that quantum mechanics does for blackbody radiation. In particular, we 

want to reduce the set of micro-states and obtain an expression for expected energy that is 

analogous to the Planck distribution. But first we must ask: What physical process acts to restrict 

the set of accessible micro-states on an earthquake fault?

The answer is that the finite strength of the fault restricts the set of micro-states. When 

we defined our basis displacement fields, we allowed the amplitudes to be arbitrarily large, 

which means that the shear stress on the fault can be arbitrarily large. But this is not possible, 

because the fault has finite strength. A micro-state cannot have shear stress that exceeds the 

fault’s yield stress, because an earthquake ould occur before such a state could be attained.

Ideally, we would require that the shear stress be less than the yield stress at every point 

on the fault, but such a requirement is not mathematically tractable. This is because the 

mathematical properties of the Fourier transform, such as Parseval’s theorem, tell us about meansquare values rather than absolute maximum values. So instead, we require that the expected 

mean-square shear stress 〈 

 

〉 be less than some fixed, finite value that depends on the strength 

of the fault.

Recall equation 25 that relates expected mean-square shear stress to expected energy:
23

⟨ total

 

⟩ 

 

 

 ∫⟨ 

 ⟩ (30)

Now we consider the value of 〈 

 

〉 to be fixed, and explore qualitatively the 

consequences for the expected energy ⟨ 

 ⟩. For small , the factor of in the integrand 

means that ⟨ 

 ⟩ contributes very little to the stress. In other words, small- states do not 

“feel” the finite strength of the fault very much. So we expect that ⟨ 

 ⟩ remains close to 

the equipartition value .

For large , the factor of in the integrand means that ⟨ 

 ⟩ contributes a lot to the 

stress. In other words, large- states do “feel” the finite strength of the fault substantially. So we 

expect that ⟨ 

 ⟩ is substantially less than the equipartition value . In fact, ⟨  ⟩

must approach zero as quickly enough so that the integral remains finite.

Finally, we note that the factor is isotropic, that is, it does not depend on the direction 

of the wave vector. So, we expect that the finite strength of the fault affects the value of 

⟨  ⟩ in an isotropic manner.

To make this explicit, we introduce an energy rolloff function which describes how 

the energy decreases with increasing wavenumber (fig. 5). It is defined by:

〈  〉 

 

 

 (31)

 (32)

 (33)

⟨ total

 

⟩ 

 

 ∫ 

 (34)

This is clearly analogous to the Planck distribution. The expected energy ⟨ 

 ⟩ is 

close to when , but approaches zero when .

Up until now, everything we did was scale-invariant, that is, there was no preferred 

length scale. But now, for the integral of 

 

to be finite, has to approach zero over 

some characteristic length, . We call this the correlation length. To make the definition precise, 

we define the correlation length to be ̂, where ̂ is chosen to maximize the value of 

 ( ̂) ̂ . (Note that this definition may not agree ith use of the term “correlation length” as it 

appears elsewhere in the literature.)

So we have shown the following: If the fault has finite strength, then there must be a 

length scale in the system, which is the correlation length.

The integral ∫ 

 must be proportional to  

, because the integrand is 

significant only when is less than some multiple of , and is of order unity. So, if we 

regard 〈 

 

〉 as being fixed by the strength of the fault, the last equation above implies a relation 

between the thermodynamic parameter and the correlation length:

 

 

 (35)
24

We see from this relation that choosing a value for is equivalent to choosing a 

correlation length.

Note that the energy rolloff function is isotropic (that is, we write and not 

 ) because we expect the finite fault strength to have an isotropic effect as mentioned 

earlier. It then follows that the expected power spectral density, 〈 

 

  〉, is also isotropic:

⟨ rms

 

  ⟩ 

 

 

〈  〉 

 

 

 (36)

Now we need to derive a formula for . Because we know that , the main 

task is to find a formula for the asymptotic behavior of as . 
25

Scaling Relations and Self-Similarity

We derive the asymptotic behavior of the energy rolloff function from a set of 

earthquake scaling relations. First, we state four scaling relations, and obtain a formula for the 

total area of stress perturbations of various sizes. Then, we can calculate the asymptotic behavior 

of by requiring it to be consistent with the area formula.

A key assumption is that earthquakes are self-similar or scale-invariant. Self-similarity is 

supported by observations such as the independence of stress drop and earthquake size, and the 

fractal geometry of faults. It is sometimes stated that if earthquakes are self-similar, then the 

stress distribution ought to be self-similar too. But a self-similar stress distribution implies a 

 

 power spectral density, which, as we have noted, leads to infinitely large stresses (Andrews 

and Barall, 2011). Here we argue that, even though earthquakes are self-similar, the stress 

distribution is not self-similar. This deviation arises from the fact that earthquakes of different 

sizes repeat at different rates.

Scaling Relations

We begin by following Hanks (1977, 1979) who states three empirical scaling relations 

and then derives a formula for the rate of occurrence of earthquakes as a function of rupture size.

Scaling relation 1 is the Gutenberg-Richter relation (Gutenberg and Richter, 1954). Let 

be an earthquake magnitude, and let be the rate of occurrence of earthquakes with magnitude 

 . Then

 (37)

where is some constant.

Scaling relation 2 is the relation between magnitude and seismic moment (Hanks and 

Kanamori, 1979). Let be seismic moment, which is defined as

 

 

(38)

where is fault slip, and is the radius of the fault rupture area. Then

 (39)

where is some constant. Notice that the implied constant of proportionality in the definition of 

 can be absorbed into the constant .

Scaling relation 3 is the statement that earthquakes are self-similar or scale-invariant. It is

 (40)

If rock behaves like a linear elastic material, then stress drop is proportional to . 

So an equivalent way to state the assumption of self-similarity is that stress drop is independent 

of earthquake size, which we write as

 

 

 (41)
26

Algebraically combining the above three scaling relations (see appendix C) shows that 

the rate of occurrence of earthquakes with radius between and is:

 

 

 

 (42)

Hanks suggests that in the distribution of shear stresses, the number of positive stress 

perturbations of radius should be proportional to , on the basis of the following 

conceptual model. “[I]magine a planar fault surface large in comparison to any earthquake 

source dimension of interest, and a population of incipient earthquakes to occur upon it; the 

earthquake population is characterized by the frequency of occurrence [ ] ... Before any of 

the earthquakes occur, all of the stress differences that will be realized at the time of occurrence 

for each and every event exist on the fault surface in ‘potential’ form” (Hanks, 1979).

There is a problem with this interpretation of . Suppose we try to calculate the 

total area of all the stress perturbations with radius :

 ∫ 

 

 

 

 ∫ 

 

 

 

 (43)

The integral diverges. The divergence tells us that there isn’t enough room on the fault surface to 

hold all the small stress perturbations if they follow Gutenberg-Richter statistics. This is yet 

another example of an ultraviolet divergence because the divergence occurs at small scales.

The problem arises because Gutenberg-Richter gives the rate at which 

earthquakes occur over an interval of time, while the stress perturbations exist simultaneously at 

an instant of time, so it’s incorrect to equate the t o.

To find the correct statistics, we must take account of the fact that earthquakes repeat, 

and smaller earthquakes repeat more frequently. Because earthquakes of different sizes repeat at 

different rates, the stress perturbations existing at an instant of time will have statistics that differ 

from Gutenberg-Richter. So we introduce one more scaling relation.

Scaling relation 4 is the statement that slip deficit accumulates at a constant rate. It is

 (44)

where is the amount of time it takes for an earthquake to repeat. This relation expresses the 

idea that all points on the fault must have the same average slip rate.

The number of earthquakes that occur over an interval of time is proportional to . So 

we can write

 

 

 

 P (45)

where is the number of stress perturbations with radius between and that exist 

at an instant of time. Combining this with the formula for and the scaling relations 

 , we get
27

 P 

 (46)

Then the total area of all stress perturbations with radius is:

 ∫ P 

 

 

 

 ∫ 

 

 

 (47)

By taking earthquake repetition into account, we have resolved the ultraviolet divergence.

Energy Rolloff and Power Spectrum

The scaling relation tells us about the size distribution of small positive 

stress perturbations. Since small stress perturbations are associated with large-wavenumber or 

short-wavelength Fourier modes, the scaling relation also tells us something about the expected 

power spectral density 〈 

 

  〉 of the stress field, in the limit of large wavenumber. So we 

ask the following question: what asymptotic behavior of the power spectral density, in the limit 

of large , is consistent with the scaling relation?

Answering that question requires a lengthy and technical calculation, which is given in 

appendix D. Here we state the result. To be consistent with the scaling relation, the asymptotic 

behavior of the expected power spectral density is

⟨ rms

 

  ⟩ 

 

 

  

as (48)

The corresponding asymptotic form of the energy rolloff function is

 

 

 

  

as (49)

The parameter has dimensions of length, and is a length scale related to the correlation length. 

The occurrence of a length scale should not be surprising, because we have previously shown 

that the finite strength of the fault implies the existence of a length scale.

It should be noted that to derive asymptotic formulas for the power spectral density and 

energy rolloff, one must make additional assumptions about the statistical properties of the stress 

field. For the formulas above, we assumed that the high-wavenumber part of the stress field has 

Gaussian statistics. The formulas can change if different assumptions are made. This is explored 

in appendix D.

The power spectral density differs from the 

 

that would be obtained by assuming the 

stress itself is self-similar, although the difference is the slowly growing factor 

 

. The 

formula resolves our ultraviolet divergence, because for any ,

∫

 

 

  

 

 

 

 

 

 

 (50)

So the high- basis displacement fields no longer contribute an infinite amount of stress.
28

State Restriction and Probability Distribution

Our plan for resolving the ultraviolet divergence is to restrict the set of accessible microstates, so that the expected mean-square stress is finite. This requires that the expected energy 

rolls off to zero as the wavenumber becomes large. In the previous section, we found the 

asymptotic behavior of the energy rolloff. In this section, we complete the plan by showing how 

to restrict the set of states and generate random values in accordance with the resulting 

probability distribution.

Energy Rolloff Function

In the last section we found the asymptotic behavior of the energy rolloff function, ,

in the limit of large , and we also know that . We need to choose a functional form for 

 that has this behavior. This could be done in any number of ways, but we have chosen the 

following form because it is among the simplest formulas with correct asymptotic behavior:

 

 

 

  

 (51)

The constants and are adjustable parameters with dimensions of length. Generally, we take 

 , and adjust to obtain the desired correlation length. Appendix E contains a recommended 

numerical algorithm. Recall that the correlation length is defined to be ̂, where ̂ is 

chosen to maximize the value of ( ̂) ̂ . 

In addition to describing the decrease in energy with increasing wavenumber, the energy 

rolloff function also establishes a relation between the thermodynamic parameter and the 

expected mean-square stress ⟨ 

 

⟩. Given one, the other can be calculated using

⟨ total

 

⟩ 

 

 ∫ 

 

 

 

(52)

Appendix E contains a recommended numerical algorithm for evaluating the improper integral.

State Restriction

For each wave vector , we impose a maximum allowed amplitude 

for the basis displacement field:

 

 max  (53)

Because energy is proportional to the square of the amplitude, the restriction can 

equivalently be described as imposing a maximum allowed energy for the basis 

displacement field:

 

 max  (54)

Subject to this restriction, the amplitude and phase are chosen 

randomly according to the Boltzmann probability distribution for a system with two degrees of 
29

freedom. Each or is chosen in such a way so that the expected energy 

satisfies 〈 〉 . It can be shown that this condition is satisfied if we define

 max 

 

 

 (55)

where is the solution to the equation

 

 

 

 (56)

Furthermore, if is a random variable uniformly distributed between 0 and 1, then

 

 

 

 (57)

is randomly distributed according to the Boltzmann distribution for a system with two degrees of 

freedom, obeys the restriction 

  , and has the expected energy 

〈 〉 . Appendix E contains proofs of these formulas, plus a recommended 

numerical algorithm for computing . This provides a practical method for calculating random 

energy values that obey the Boltzmann distribution, our state restriction, and our energy rolloff.
30

Algorithm to Generate the Stress Field

In this section we present the algorithm for generating a random stress field. It is the sum 

of a set of plane waves, each associated with a basis displacement field. Amplitudes and phases 

are randomly selected, independently for each wave vector, according to the Boltzmann 

probability distribution, with the amplitudes restricted to conform with the energy rolloff 

function.

The stress field produced by this algorithm is a “ra ” field, hich generally must be 

post-processed before being used in a dynamic rupture simulation. Later sections will discuss 

post-processing operations, which typically include the following: limiting the extent of the 

simulated rupture; scaling the stress field to conform to the assumed fault friction; and filtering 

the stress field to suppress high-frequency components that cannot be resolved by the dynamic 

rupture simulation code. For models where the initial normal stress varies over the fault surface, 

the “ra ” values must be scaled proportionately to the normal stress; thus, the “ra ” values 

should be regarded as being in “friction coefficient space.”

The algorithm consists of three parts. First, there are several preliminary steps. Second, in 

the main part of the algorithm, a set of random amplitudes and phases are generated. Third, the 

corresponding stress fields are summed to produce the final result.

Preliminary Steps

Step 1. Choose a correlation length, . The next section will provide guidance in 

choosing the correlation length. The value of is typically about equal to the diameter of the 

earthquake rupture.

Step 2. Solve for the parameters and that appear in the energy rolloff function , 

as described in the previous section. Typically we choose .

Step 3. Choose a value for the thermodynamic parameter . If the stress field is going to 

be scaled during post-processing, then may be chosen arbitrarily because it plays the role of a 

scale factor which is adjusted during post-processing. Otherwise, could be chosen by assuming 

a value for the mean-square stress and using the relation between and mean-square stress given 

in the previous section.

 tep 4. Choose the size of the “ra ” fault surface, hich is vertically and alongstrike. For reasons explained in the next section, the “ra ” fault surface is typically much larger 

than the final fault surface used in a dynamic rupture simulation. Typically we choose .

Step 5. Lay out a grid of points on the fault surface, which are the locations where the 

stress field is computed. Let and be the number of grid points vertically and along-strike, 

respectively. We require that and be powers of (so that we can use the simplest form of 

the Fast Fourier Transform), and typically we choose . The grid points are

 

 for (58)

 

 for 

(59)

where and are integers. The grid is considered to satisfy periodic boundary conditions. 

Corresponding to this grid of points, there is a grid of wave vectors which is

 for (60)
31

 for (61)

where and are integers. The Nyquist wavenumbers are defined to be the highest 

wavenumbers in the grid,

 

 

(62)

 

 (63)

 tep 6. elect a set of “active” ave vectors as follows. For each pair of wave 

vectors and that satisfy the condition

 

 

 

 

 

 

 (64)

select one member of the pair to be “active” (it doesn’t matter hich one is selected). We only 

include one member of each pair because and correspond to the same 

physical state. Notice that in the typical case where 

  

, the above condition describes a 

circle in the -plane. We use a circle rather than a square so that the total stress field will be 

isotropic; a square allows shorter wavelengths along the grid diagonals than along the grid axes. 

We exclude the zero wave vector because it corresponds to uniform fault slip, which produces no 

stress. We exclude the wave vectors 

 

 , 

 

 , and 

 

 

 

 because they each introduce 

only one degree of freedom, whereas our analysis assumes two degrees of freedom per wave 

vector.

Generating Random Amplitudes and Phases

The follo ing steps are performed once for each “active” ave vector .

Step 7. Calculate the energy rolloff function .

Step 8. Numerically solve the following equation for the energy truncation value :

 

 

 

 (65)

Step 9. Generate two independent random numbers and , each of which is uniformly 

distributed between and . Calculate the energy and phase according to

 

 

 

 

 (66)

 

 (67)

Note that it is not necessary to calculate the displacement field amplitude because the 

stress field can be computed directly from the energy and phase.
32

Summing the Random Stress Fields

Step 10. For each grid point , calculate the total stress field as

 total 

 ∑ √

 

 

 

 ( 

 )

 

(68)

 here the sum runs over all “active” ave vectors . The summation can be performed 

efficiently using the Fast Fourier Transform. This is our final “ra ” stress field.

Note that the above summation formula will be modified slightly in the next section, as 

part of controlling the rupture extent.
33

Controlling Rupture Extent

When performing a dynamic rupture simulation, some method must be used to limit the 

extent of the rupture. Details of several methods can be found in the SCEC/U.S. Geological 

Survey dynamic rupture code verification project (Harris, Barall, Archuleta, and others, 2009; 

Harris, Barall, Andrews, and others, 2011). The simplest technique is to impose a strength 

barrier, beyond which the fault is not permitted to rupture. Aside from being an arbitrary limit, 

the strength barrier also brings the rupture to a “hard stop,” hich can generate strong 

reflections. Another way to limit the extent of the rupture is to modify the initial stress or the 

friction properties outside the desired rupture area to make rupture propagation unfavorable. If

the modification is gradual, it can bring the rupture to a “soft stop.” For example, Andre s and 

Barall (2011) show how to create low initial stresses outside the desired rupture area by 

modifying several of the Fourier coefficients of the stress field.

Here e introduce a ne method for limiting rupture extent, hich creates a “soft stop” 

by having low initial stresses outside the rupture area. Unlike previous methods, the change in 

stress arises spontaneously from the randomly generated heterogeneous initial stress field and is 

neither arbitrarily imposed nor created by manipulation. The idea is to generate a random stress 

field over a very large area, much larger than the actual fault. Then, the large area is searched to 

find a fault-sized subarea which, just by chance, has high stress in the center and low stress at the 

sides.

Figure 6 shows examples of stress distributions produced by our method. Rupture extent 

is limited by areas of low initial stress located left, right, and below the rupture area. High initial 

stress prevails within the rupture area. The shape of the rupture area is random, but its size is 

approximately equal to the correlation length. We will see that this provides a method for 

selecting the correlation length.

Depth Conditioning

Our focus here is on earthquakes large enough to rupture the full depth of the 

seismogenic zone. The base of the seismogenic zone is defined by a systematic change in the 

properties of the fault, such as a transition from brittle to ductile behavior. It is therefore 

appropriate to introduce a non-random method to limit rupture extent at the bottom of the fault.

Rupture extent along-dip is controlled by depth conditioning, a post-processing step that 

systematically reduces the shear stress with increasing depth. The shear stress is multiplied by a 

depth-conditioning function, , which varies as a function of depth:

 

   (69)

This step is performed after the stress has been scaled and filtered. (Scaling and filtering 

are described later.) In the examples shown here, the depth conditioning function is taken to be 

piecewise linear:

 

 {

 for m

 m for m m

 m for m m

(70)
34

It is, of course, possible to make the depth conditioning function be smooth rather than 

piecewise linear.

Selection Process

Rupture extent along-strike is controlled by a selection process, which selects the fault 

surface from a randomly generated stress distribution thousands of times larger.

The selection process is illustrated in figure 7. In this example, the fault surface is 48.0 

km along-strike and 19.5 km along-dip, and grid points are spaced 75 m apart. We generate the 

stress distribution using , which implies that the large stress distribution is a 

square measuring 2,457.6 km on a side. With these choices, the area of the large stress 

distribution is more than 6,000 times larger than the area of the fault. We consider every possible 

position of the fault within the large stress distribution, and pick the position which is optimal in 

the sense of having high stress in the center of the fault and low stress at the sides of the fault. 

We then cut out the portion of the stress distribution at the optimal position, and use it as the 

stress distribution over the fault surface.

Remark: The value of and is determined mainly by the available computer 

memory. The value 32,768 allows the algorithm to be performed on a computer with 48 

gigabytes of memory.

Two improvements can be made to the selection process. We described the selection 

process as looking for high stress in the center of the fault and low stress at the sides, but the 

reverse would be just as good. If we found a location with low stress in the center of the fault and 

high stress at the sides, we could simply change the sign of the stress, because our method of 

randomly generating the stress field is symmetric with respect to a change of sign. This 

improvement effectively doubles the area of the large stress distribution, for no extra work.

A further improvement is to allow the stress field to be complex-valued. Then, instead of 

allowing just a change of sign, we can allow an arbitrary change in the complex phase. (Note that 

changing the complex phase by 180 degrees is equivalent to a change in sign.) This works 

because our method of generating the stress field is symmetric with respect to a change in the 

complex phase, and the result is another effective doubling of the area of the large stress 

distribution. Our algorithm incorporates this improvement.

Selection Algorithm

We now present an algorithm for finding the location of the fault within the large stress 

distribution. A mathematical justification of the algorithm is given in appendix F.

Step 1. Create a selection template. It is a function 

 which is positive in the 

portion of the fault surface where we hope to have high stress and negative in the portion of the 

fault surface here e hope to have lo stress. It is zero in areas here e don’t care about the 

stress. Our objective is to find the place in the large randomly generated stress distribution that is 

most similar to the selection template.

Note that is defined for all points throughout the large stress 

distribution, including points that are outside the fault surface. It is created under the assumption 

that the upper left corner of the fault surface is at the origin 

 . Generally 

 would be zero outside the fault surface, although nonzero values may extend a short 

distance outside the fault surface to allow the selection template to taper to zero rather than have 

a discontinuity at the edge of the fault surface. To make the selection template sensitive to 
35

variations in stress rather than the absolute level of stress, it is adjusted to have zero DC 

component, that is,

∑ 

 

 (71)

where the sum runs over all points in the large stress distribution. An example of a 

selection template is shown in figure 8. In the example, is positive in the center 30.0 

km of the fault surface, because we aim to produce a rupture about 30-km wide. It is negative on 

the left and right sides of the fault surface. The selection template is zero in the bottom 2.5 km of 

the fault surface, because depth conditioning forces the stress to be lo there, so e don’t care 

about the randomly generated stresses near the bottom of the fault surface.

Step 2. Calculate a complex-valued random shear stress field:

 total

 

  ∑ √

 

 

 

 ( 

 ) 

 

(72)

The randomly generated energies and phases are the same ones described in 

the previous section. In fact, the previously described is just the real part of 

 

 

 . The summation can be done efficiently using the Fast Fourier Transform.

Step 3. Compute the convolution of the selection template and the complex-valued stress 

field:

 total

 

  ∑  total

 

 

 

 

(73)

The convolution can be computed efficiently using the Fast Fourier Transform. Because we use 

periodic boundary conditions on the large stress distribution, the sums and on 

the right-hand side are taken to wrap around cyclically. Note that convolution is usually defined 

using and on the right-hand side, but for our purposes it is more convenient to 

use the plus sign.

Step 4. Find the point ̃ ̃ in the large stress distribution that maximizes the absolute 

value of the convolution, | 

 

 ̃ ̃  |. This will be the upper left corner of the fault 

surface within the large stress distribution.

Step 5. Let be the complex phase of the convolution 

 

 ̃ ̃  at the point 

where its absolute value is maximized. In other words, write

 total

 

 ̃ ̃  | total

 

 ̃ ̃  | 

 

 (74)

Then calculate the real-valued random shear stress as:

 

 ∑ √

 

 

 

 

 

 

(75)
36

Note that each cosine term is phase-shifted by the same amount .

Step 6. Place the upper left corner of the fault surface at location ̃ ̃ . Cut out the part 

of the stress distribution that lies within the fault surface, and discard the rest. The 

remaining is the randomly generated stress field on the fault surface, subject to further 

post-processing.

Correlation Length

The selection process gives us a way to determine the correlation length , although it 

requires some trial and error. If the correlation length is set too small, then it will be impossible 

to find a stress pattern that supplies high stress over the desired area of the rupture. As a result, it 

 on’t be possible to successfully execute a dynamic rupture simulation.

Our procedure is to select the smallest correlation length that reliably gives stress patterns 

that can sustain a rupture of the desired size. Empirically, it appears that the correlation length 

should be about equal to the desired diameter of the rupture.

Figure 9 shows the effect of varying the correlation length. It shows stress patterns 

generated with correlation lengths of 10, 30, and 100 km, with a desired rupture size of 30 km. A 

correlation length of 10 km is too small, because it is impossible to find any high-stress area 

large enough to sustain a rupture 30 km across. A correlation length of 100 km is larger than 

necessary and yields a fairly clean division of the fault into high and low stress areas. By 

comparison, a correlation length of 30 km is about right. It yields a high stress area large enough 

to sustain a rupture 30 km across, yet the edges of the high-stress area are ragged, with the high 

and low stress areas intruding into each other.

To carry out this procedure, we suggest generating stress distributions with several 

different correlation lengths, beginning with the desired rupture diameter, and then adding 

correlation lengths both larger and smaller by factors of two or three. Then, check which 

correlation lengths produce stress patterns able to sustain ruptures of the desired size. Correlation 

lengths that are too small can often be identified by simple visual inspection of the stress pattern, 

noting the absence of any high-stress area approximating the desired rupture size, as shown in 

the 10-km panel of figure 9. It is also possible to identify too-small correlation lengths by 

attempting to perform a dynamic rupture simulation on the stress pattern and noting that the 

rupture either cannot be nucleated, or dies out well before it reaches the desired size. The 

smallest correlation length that reliably leads to successful dynamic rupture simulations of the 

desired size is the correlation length that should be used.

Application to Smaller Earthquakes

A notable feature of the method just described is that the rupture is limited vertically and 

horizontally by two different techniques. The rupture is limited at the bottom by imposing depth 

conditioning, whereas the rupture is limited at the left and right by a selection process.

The reason for this difference is that our discussion focuses on the case where the 

earthquake is large enough to rupture the full depth of the seismogenic zone. In nature, large 

earthquake ruptures stop at the bottom of the seismogenic zone because of a change in the fault 

properties, specifically, a transition from brittle to ductile behavior. What causes such ruptures to 

stop as they propagate along-strike is not well understood, but does not appear to be the result of 

any major change in fault properties. Our use of different techniques in the along-dip and alongstrike directions reflects this difference in nature.
37

In the case of smaller earthquakes, where the rupture does not extend all the way from the 

top to the bottom of the fault, it would not be appropriate to use two different techniques. In this 

case, depth conditioning should not be used. Instead, the selection process alone is used to limit 

rupture extent in all directions. The selection template is chosen to be positive in a 

region that approximates the desired rupture extent, which could be a circle, square, ellipse, or 

rectangle.
38

Post-Processing

The previous section completed the construction of a raw random stress field on a fault 

surface. Before the stress field can be used as input to a dynamic rupture simulation, it is 

generally necessary to perform some post-processing steps. The steps required will depend on 

the details of the model, such as the friction parameters, the variation of normal stress with 

depth, and the characteristics of the simulation software. A typical post-processing sequence 

could include the following steps:

 Stress field rotation. If stresses are higher near the Earth’s surface than in the lo er 

part of the fault, then rotate the stress pattern 180 degrees to place the higher stresses 

in the lower part of the fault. This is useful because high stress at depth is usually 

more effective at driving the rupture.

 Filtering. Apply a low-pass spatial filter to suppress stress fluctuations that are below 

the resolution of the dynamic rupture simulation code.

 Scaling. Apply a scale factor and offset (that is, a linear transformation) to the stress 

pattern. Scaling is used to make the initial stresses consistent with the friction law. In 

the high-stress (center) portion of the fault surface, the stresses must be high enough 

to sustain the rupture but belo the fault’s yield stress. In the lo -stress (side and 

bottom) portion of the fault surface, the stresses must be low enough to stop the 

rupture.

 Hypocenter selection. Select a point within a region of relatively high initial stress to 

be the hypocenter.

 Depth conditioning. As discussed earlier, stresses at the bottom of the fault are 

reduced, to limit the extent of the rupture along-dip.

 Scaling by normal stress. The raw stress values are multiplied by a factor proportional 

to the initial normal stress on the fault. This is done because, in general, friction 

forces are proportional to the normal stress. Another way of saying this is that we 

treat our randomly generated values as if they lie in “friction coefficient space.”

Appendix G contains algorithms and further details for post-processing.
39

Nucleation

A dynamic rupture simulation must have some way to start, or nucleate, the rupture. If 

the initial conditions are set up so that the initial shear stress is every here less than the fault’s 

yield stress, then nothing will happen unless some artificial means of nucleating the rupture is 

implemented. Descriptions of several nucleation methods can be found, e.g., in the SCEC 

Dynamic Rupture Code Validation benchmarks (Harris, Barall, Archuleta, and others, 2009; 

Harris, Barall, Andrews, and others, 2011), and in Bizzarri (2010).

Here we introduce a new two-stage nucleation method. The method is designed to let the 

rupture respond to the heterogeneous initial stresses beginning very early in the simulation. The 

first stage is a small zone of forced rupture surrounding the hypocenter, which starts the rupture 

going. The second stage is a larger zone of reduced fracture energy, which allows the rupture to 

become large enough so that further expansion is possible without additional assistance.

The Problem of Nucleation

One of the challenges in performing dynamic rupture simulations is to find a way to 

nucleate the rupture. The problem is that, for the range of parameter values typically used in 

dynamic rupture simulations, the rupture does not become self-sustaining until it has a radius of 

several kilometers (this is discussed further below). So, one must introduce an artificial 

mechanism to propel the rupture forward until it becomes large enough to sustain itself.

We do not attempt to give a comprehensive treatment of nucleation techniques that have 

been used in the past, but refer the reader to the SCEC website (http://scecdata.usc.edu/cvws) for 

detailed descriptions of several techniques.

One simple technique is to establish a nucleation zone surrounding the hypocenter, with a 

radius of several kilometers, and to set the shear stress within the nucleation zone to be higher 

than the fault’s yield stress. This high stress can be applied instantaneously at the beginning of 

the simulation or gradually over some interval of time; both variants have been used in the SCEC 

benchmarks. This technique has the disadvantage that the final slip in the nucleation zone is 

much larger than elsewhere on the fault. Also, it injects a large amount of energy into the 

rupture, which can make earthquake characteristics sensitive to the details of the nucleation 

process.

Another technique is to alter the frictional properties of the fault near the hypocenter, for 

example, by lowering the yield stress to below the initial shear stress. Again, this can be done 

either instantaneously at the start of the simulation or gradually over an interval of time. Often, 

the altered frictional properties are applied in an expanding circle surrounding the hypocenter, 

creating a forced rupture front (Andrews, 2004). This technique has the disadvantage that rupture 

propagation is controlled by the pre-determined nucleation process, at least until the rupture is 

large enough to become self-sustaining. Note that by the time the rupture becomes selfsustaining, it already fills an appreciable fraction of the fault surface.

We are interested in studying the effects of stochastically generated heterogeneous initial 

stresses. So, our interest is to have the rupture propagation and the slip pattern be controlled by 

the heterogeneous initial stresses, and not by the nucleation process, to the greatest extent 

possible. The new nucleation technique described here is an attempt to reduce the influence of 

the nucleation process on the rupture propagation and slip pattern, as compared to prior 

techniques.
40

The Day Radius

Day (1982) gives a formula for the minimum radius that a rupture must attain, so that 

further expansion of the rupture is energetically favorable, that is, so that the elastic energy 

released by expansion of the rupture exceeds the fracture energy consumed:

 

 

 

 

 

 

 (76)

Here, is the shear modulus, is the initial shear stress, is the static yield stress, is the 

dynamic sliding stress, and is the slip-weakening critical distance, that is, the amount of slip 

over which the frictional stress declines from to . The fracture energy per unit area is 

  .

The above formula for the Day radius was derived for the case of uniform shear stress, 

for a circular rupture zone, on a fault with uniform frictional properties, embedded in a uniform 

linear elastic fullspace. However, we have found the formula to be a useful approximation in 

other cases. When stresses or friction parameters vary with position on the fault, we insert into 

the Day formula their average values within a neighborhood of the hypocenter.

For the range of parameter values typically used in dynamic rupture simulations, the Day 

radius is on the order of 3 km. The nucleation problem is that a simulated rupture does not 

become self-sustaining until its radius is approximately the Day radius.

Notice that is proportional to the slip-weakening critical distance . If is made to 

vary with location on the fault surface, then and the fracture energy are also variable. In 

particular, if is made proportional to the distance from the hypocenter, then it is possible for 

the rupture radius to be approximately equal to over an extended range of rupture sizes. Note 

that only the value of near the rupture front is relevant to the question of whether or not it is 

energetically favorable for the rupture to expand, because it is only near the rupture front that 

fracture energy is consumed. These observations are the basis of our two-stage nucleation 

method.

(In principle, the slip-weakening critical distance is a physical parameter. In practice, 

the value of used in a dynamic rupture simulation is much larger than the critical distance 

observed in laboratory experiments. It is chosen large enough so that the model is able to resolve 

the cohesive zone, which is the area where slip has begun but not yet reached the critical 

distance. It is for this reason that we are free to manipulate to help achieve nucleation. Also, 

notice that as computers become more powerful and are able to run models with finer resolution, 

the practical value of the Day radius will decrease.)

Two-Stage Nucleation Method

Our two-stage method constructs two circular zones surrounding the hypocenter, a small 

zone of forced rupture with radius , and a larger zone of reduced fracture energy with radius 

 . We take equal to one-fourth of . The radius is approximately equal to the 

Day radius, although it may vary somewhat depending on the model parameters, and some 

experimentation may be necessary to choose a value for that produces reliable nucleation. 

For example, in the SCEC Dynamic Rupture Code Validation benchmarks TPV18 through 

TPV21, the radii are =900 m and =3,600 m. Appendix G contains additional 

information on selecting and .
41

The first stage occurs within the zone of forced rupture. Each point within the zone is 

assigned a time, , of forced rupture,

 {

 

 if F 

 F  F  if F F 

(77)

where is distance from the hypocenter, and is the shear-wave velocity. At the time , the 

frictional stress drops immediately to the dynamic sliding stress. The result is a circular front of 

forced rupture, which expands outward from the hypocenter at a speed of when 

 , and when . Forced rupture does not occur outside the 

zone.

The first stage by itself is not sufficient to nucleate the earthquake, because the zone of 

forced rupture is much smaller than the Day radius. So we introduce a second stage within the 

zone of reduced fracture energy. Each point within the zone is assigned a reduced value 

 of the 

slip-weakening critical distance,

 

 {

  if F 

 F 

 if F F 

(78)

where is the value of the slip-weakening critical distance outside the zone of reduced fracture 

energy. Notice that 

 

is proportional to when . As discussed above, the 

result is a variable Day radius, which permits the rupture radius to be approximately equal to the 

Day radius over the size range . In other words, the radius at which 

rupture expansion becomes energetically favorable is reduced by approximately a factor of 10.

Our two-stage nucleation method has been used successfully in six of the SCEC Dynamic 

Rupture Code Validation benchmarks.
42

Conclusions

This report presents a new method for generating random heterogeneous initial stresses, 

for use in dynamic rupture simulations. The method employs concepts and techniques drawn 

from thermodynamics and statistical mechanics.

A pattern of slip on the fault surface is taken to be analogous to a micro-state of a 

thermodynamic system. Fault slip distorts the surrounding elastic medium, and the resulting 

elastic energy is taken to be the energy of the micro-state. Then, the Boltzmann probability 

distribution gives the probability of each possible pattern of slip and corresponding pattern of 

stress.

Although the idea is simple in principle, there are challenges in carrying it out. To make 

the method computationally feasible, the system must be decomposed into a set of degrees of 

freedom, so that the Boltzmann distribution can be applied independently to each degree of 

freedom. We decompose the system by introducing a set of three-dimensional basis displacement 

fields, each of which is an approximate solution to the elastostatic problem for a sinusoidal 

pattern of fault slip. The basis fields have the crucial property of orthogonality, which means that 

the energy of a linear combination of basis fields equals the sum of the energies of the individual 

fields. Orthogonality implies that the basis displacement fields behave like the degrees of 

freedom of a thermodynamic system, allowing the amplitudes of the basis fields to be selected 

independently according to the Boltzmann distribution.

A straightforward application of the Boltzmann distribution leads to infinitely large 

stresses on the fault. This phenomenon is a consequence of the equipartition theorem and is 

called an ultraviolet divergence because the divergence occurs in the limit of small wavelength.

The finite strength of the fault acts to restrict the possible states of the system, because a 

state that includes shear stresses in excess of the fault’s yield stress cannot occur. When the finite 

fault strength is taken into account, the energy spectrum must deviate from equipartition, with 

the expected energy falling off to zero as the wavenumber approaches infinity. To find a formula 

for the energy falloff, we reconsider several earthquake scaling relations—the Gutenberg-Richter 

relation, the magnitude-moment relation, and earthquake self-similarity—and we introduce a 

new scaling relation of constant slip deficit accumulation. From these scaling relations, we are 

able to deduce that the expected power spectral density in the limit of large wavenumber is 

proportional to 

 

 

 

, where is the scalar wavenumber and is a length scale related 

to the correlation length. With this formula in hand, the ultraviolet divergence is resolved, and 

we are able to construct a complete algorithm for generating random patterns of stress. It is 

noteworthy that the algorithm has only one free parameter, which is the correlation length.

We have also introduced a new selection method for limiting the extent of a simulated 

earthquake rupture. We generate a random stress pattern over an area thousands of times larger 

than the fault surface and then select a fault-sized portion of the stress pattern which, just by 

chance, has a positive stress perturbation of the desired size. The simulated rupture then stops 

automatically when it reaches the low-stress areas at the edge of the stress perturbation, without 

the need to impose hard barriers.

In addition to limiting rupture extent, the selection method also yields one-point statistics 

for the stress distribution that are heavy-tailed (by comparison to Gaussian). Furthermore, the 

selection method provides a way to determine the appropriate correlation length, because if the 

correlation length is made too small then it is impossible to find a positive stress perturbation of 

the desired size. With this observation, the number of free parameters governing the probability 

distribution is reduced to zero. (Of course, there remain many free parameters in the post-
43

processing procedure, nucleation process, friction law, and elsewhere in the dynamic rupture 

simulation.)

Finally, we have introduced a new two-stage nucleation method, which combines a small 

zone of forced rupture with a larger zone of reduced fracture energy. The two-stage nucleation 

method is designed so that the rupture propagation can become responsive to the heterogeneous 

initial stresses beginning very early in the simulation.

All of this has been successfully tested in the Southern California Earthquake Center 

Dynamic Rupture Code Validation group. Ten different dynamic rupture codes performed 

simulations using heterogeneous initial stress patterns generated by our method. All ten codes 

produced results in excellent agreement with each other, demonstrating that our method is wellsuited to generating heterogeneous initial stress conditions for dynamic rupture simulations.

Acknowledgments

This work is supported by the U.S. Geological Survey Earthquake Science Center, by the 

Southern California Earthquake Center, and by the Pacific Gas and Electric Company. Thanks to 

Joe Andrews and Brad Aagaard for reviewing the manuscript.
44

References

Aagaard, B.T., and Heaton, T.H., 2008, Constraining fault constitutive behavior with slip and 

stress heterogeneity: Journal of Geophysical Research, v. 113, B04301, 18 p., 

doi:10.1029/2006JB004793.

Abaimov, S.G., 2009, Critical behavior of slider-block model: arXiv:0902.3767v4 [condmat.stat-mech], available at http://arxiv.org/abs/0902.3767.

Andrews, D.J., 1974, Evaluation of static stress on a fault plane from a Green's function: Bulletin 

of the Seismological Society of America, v. 64, no. 6, p. 1,629–1,633.

Andrews, D.J., 1980, A stochastic fault model, 1. Static case: Journal of Geophysical Research, 

v. 85, no. B7, p. 3,867–3,877, doi: 10.1029/JB085iB07p03867.

Andrews, D.J., 1981, A stochastic fault model, 2. Time-dependent case: Journal of Geophysical 

Research, v. 86, no. B11, p. 10,821–10,834, doi: 10.1029/JB086iB11p10821.

Andrews, D.J., 2004, Rupture models with dynamically determined breakdown displacement: 

Bulletin of the Seismological Society of America, v. 94, no. 3, p. 769–775, doi: 

10.1785/0120030142.

Andrews, D.J., and Barall, M., 2011, Specifying initial stress for dynamic heterogeneous 

earthquake source models: Bulletin of the Seismological Society of America, v. 101, no. 5, p. 

2,408–2,417, doi: 10.1785/0120110012.

Beroza, G.C., and Mikumo, T., 1996, Short slip duration in dynamic rupture in the presence of 

heterogeneous fault properties: Journal of Geophysical Research, v. 101, no. B10, p. 22,449–

22,460, doi: 10.1029/96JB02291.

Bizzarri, A., 2010, How to promote earthquake ruptures—Different nucleation strategies in a 

dynamic model with slip-weakening friction: Bulletin of the Seismological Society of 

America, v. 100, no. 3, p. 923–940, doi: 10.1785/0120090179.

Bizzarri, A., 2011, On the deterministic description of earthquakes: Reviews of Geophysics, v.

49, RG3002, 32 p., doi: 10.1029/2011RG000356.

Boatwright, J., 1981, Quasi-dynamic models of simple earthquakes—Application to an 

aftershock of the 1975 Oroville, California, earthquake: Bulletin of the Seismological Society 

of America, v. 71, no. 1, p. 69–94.

Boltzmann, L., 1876, Über die natur der gasmoleküle (On the nature of gas molecules): Wiener 

Berichte, v. 74, p. 553–560.

Bouchon, M., 1997, The state of stress on some faults of the San Andreas system as inferred 

from near-field strong motion data: Journal of Geophysical Research, v. 102, no. B6, p. 

11,731–11,744, doi: 10.1029/97JB00623.

Das, S., and Aki, K., 1977, Fault plane with barriers—A versatile earthquake model: Journal of

Geophysical Research, v. 82, no. 36, p. 5,658–5,670, doi: 10.1029/JB082i036p05658.

Day, S.M., 1982, Three-dimensional simulation of spontaneous rupture—The effect of 

nonuniform prestress: Bulletin of the Seismological Society of America, v. 72, no. 6, p. 1,881–

1,902.

Day, S.M., Yu, G., and Wald, D., 1998, Dynamic stress changes during earthquake rupture: 

Bulletin of the Seismological Society of America, v. 88, no. 2, p. 512–522.

Fefferman, C., 1971, On the convergence of multiple Fourier series: Bulletin of the 

Seismological Society of America, v. 77, no. 5, p. 744–745.
45

Guatteri, M., Mai, P.M., and Beroza, G.C., 2004, A pseudo-dynamic approximation to dynamic 

rupture models for strong ground motion prediction: Bulletin of the Seismological Society of 

America, v. 94, no. 6, p. 2,051–2,063, doi: 10.1785/0120040037.

Guatteri, M., P. Mai, P.M., Beroza, G.C., and Boatwright, J., 2003, Strong ground-motion 

predictions from stochastic-dynamic source models: Bulletin of the Seismological Society of 

America, v. 93, no. 1, p. 301–313, doi: 10.1785/0120020006.

Gutenberg, B., and Richter, C., 1954, Seismicity of the Earth and associated phenomena (2nd 

ed.): Princeton University Press.

Hanks, T.C., 1977, Earthquake stress drops, ambient tectonic stresses, and stresses that drive 

plate motions: Pure and Applied Geophysics, v. 115, no. 1–2, p. 441–458, doi: 

10.1007/BF01637120.

Hanks, T.C., 1979, Values and 

 

seismic source models—Implications for tectonic stress 

variations along active crustal fault zones and the estimation of high-frequency strong ground 

motion: Journal of Geophysical Research, v. 84, no. B5, p. 2,235–2,242, doi: 

10.1029/JB084iB05p02235.

Hanks, T.C., and Kanamori, H., 1979, A moment magnitude scale: Journal of Geophysical 

Research, v. 84, no. B5, p. 2,348–2,350, doi: 10.1029/JB084iB05p02348.

Harris, R.A., 2004, Numerical simulations of large earthquakes—Dynamic rupture propagation 

on heterogeneous faults: Pure and Applied Geophysics, v. 161, no.11/12, p. 2,171–2,181, doi: 

10.1007/s00024-004-2556-8.

Harris, R.A., Barall, M., Andrews, D.J., Duan, B., Ma, S., Dunham, E.M., Gabriel, A.-A.,

Kaneko, Y., Kase, Y., Aagaard, B.T., Oglesby, D.D., Ampuero, J.-P., Hanks, T.C., and 

Abrahamson, N., 2011, Verifying a computational method for predicting extreme ground 

motion: Seismological Research Letters, v. 82, no. 5, p. 638–644, doi: 10.1785/gssrl.82.5.638.

Harris, R.A., Barall, M., Archuleta, R., Dunham, E., Aagaard, B., Ampuero, J.-P., Bhat, H.,

Cruz-Atienza, V., Dalguer, L., Dawson, P., Day, S., Duan, B., Ely, G., Kaneko, Y., Kase, Y.,

Lapusta, N., Liu, Y., Ma, S., Oglesby, D., Olsen, K., Pitarka, A., Song, S., and Templeton, E.,

2009, The SCEC/USGS Dynamic Earthquake Rupture Code Verification Exercise:

Seismological Research Letters, v. 80, no 1, p. 119–126, doi: 10.1785/gssrl.80.1.119.

Huang. Y., Saleur, H., Sammis, C., and Sornette, D., 1998, Precursors, aftershocks, criticality 

and self-organized criticality: Europhysics Letters, v. 41, no 1, p. 43–48, doi:

10.1209/epl/i1998-00113-x.

Ising, E. (1925), Beitrag zur theorie des ferromagnetismus (Contribution to the theory of 

ferromagnetism): Zeitschrift für Physik, v. 31, no 1, p. 253–258, doi: 10.1007/BF02980577.

Kaneko, Y., Avouac, J.-P., and Lapusta, N., 2010, Towards inferring earthquake patterns from 

geodetic observations of interseismic coupling: Nature Geoscience, v. 3, p. 363–369, doi: 

10.1038/ngeo843.

Lapusta, N., 2009, Seismology—The roller coaster of fault friction: Nature Geoscience, v. 2, p.

676–677, doi: 10.1038/ngeo645.

Lavallée, D., Liu, P., and Archuleta, R.J., 2006, Stochastic model of heterogeneity in earthquake 

slip spatial distributions: Geophysical Journal International, v. 165, no. 2, p. 622–640 doi: 

10.1111/j.1365-246X.2006.02943.x.

Mai, P.M., and Beroza, G.C., 2002, A spatial random field model to characterize complexity in 

earthquake slip: Journal of Geophysical Research, v. 107, 2308, 21 p., doi: 

10.1029/2001JB000588.
46

Mavroeidis, G.P., and Papageorgiou, A.S., 2010, Effect of fault rupture characteristics on nearfault strong ground motions effect of fault rupture: Bulletin of the Seismological Society of 

America, v. 100 no. 1, p. 37–58, doi: 10.1785/0120090018.

McGarr, A., Spottiswoode, S.M., Gay, N.C., and Ortlepp, W.D., 1979, Observations relevant to 

seismic driving stress, stress drop, and efficiency: Journal of Geophysical Research, v. 84, no. 

B5, p. 2,251–2,261, doi: 10.1029/JB084iB05p02251.

Mikumo, T., and Miyatake, T., 1978, Dynamical rupture process on a three-dimensional fault 

with non-uniform frictions and near-field seismic waves: Geophysical Journal of the Royal 

Astronomical Society, v. 54, no. 2, p. 417–438, doi: 10.1111/j.1365-246X.1978.tb04267.x.

National Institute of Standards and Technology, 2011, Digital library of mathematical functions 

(2011–08–29): National Institute of Standards and Technology, http://dlmf.nist.gov.

Oglesby, D.D., and Day, S.M., 2002, Stochastic fault stress—Implications for fault dynamics 

and ground motion: Bulletin of the Seismological Society of America, v. 92, no. 8, p. 3,006–

3,021, doi: 10.1785/0120010249.

Onsager, L., 1944, Crystal statistics. I. A two-dimensional model with an order-disorder 

transition: Physical Review, v. 65, no. 3–4, p. 117–149, doi:10.1103/PhysRev.65.117.

Peyrat, S., Olsen, K., and Madariaga, R., 2001, Dynamic modeling of the 1992 Landers 

earthquake: Journal of Geophysical Research, v. 106, no. B11, p. 26,467–26,482, doi:

10.1029/2001JB000205.

Reif, F., 1965, Fundamentals of statistical and thermal physics: McGraw-Hill.

Ripperger, J., Ampuero, J.-P., Mai, P.M., and Giardini, D., 2007, Earthquake source 

characteristics from dynamic rupture with constrained stochastic fault stress: Journal of 

Geophysical Research, v. 112, B04311, 17 p., doi: 10.1029/2006JB004515.

Ripperger, J., Mai, P.M., and Ampuero, J.-P., 2008, Variability of near-field ground motion from 

dynamic earthquake rupture simulations: Bulletin of the Seismological Society of America, v.

98, no. 3, p. 1,207–1,228, doi: 10.1785/0120070076.

Royden, H., 1968, Real analysis (2nd ed.): Macmillan.

Rundle, J., Turcotte, D., Shcherbakov, R., Klein, W., and Sammis, C., 2003, Statistical physics 

approach to understanding the multiscale dynamics of earthquake fault systems, Reviews of 

Geophysics, v. 41, 1019, 30 p., doi: 10.1029/2003RG000135.

Schmedes, J., Archuleta, R.J., and Lavallée, D., 2010, Correlation of earthquake source 

parameters inferred from dynamic rupture simulations: Journal of Geophysical Research, v.

115, B03304, 12 p., doi: 10.1029/2009JB006689.

Shi, Z., and Day, S.M., 2011, 3D Simulations of dynamic rupture on rough faults: American 

Geophysical Union Fall Meeting, San Francisco, December 5–9, 2011, abstract S54C–06.

Smalley, R., Jr., Turcotte, D., and Solla, S., 1985, A renormalization group approach to the stickslip behavior of faults: Journal of Geophysical Research, v. 90, no. B2, p. 1,894–1,900, doi: 

10.1029/JB090iB02p01894.

Tullis, T.E., 2007, Friction of rock at earthquake slip rates, in Kanamori, H., ed., Treatise on

geophysics: Amsterdam and Boston, Elsevier, v. 4, chap. 5, p. 131–152.

Yamashita, T., 1976, On the dynamical process of fault motion in the presence of friction and 

inhomogeneous initial stress—Part I. Rupture propagation: Journal of Physical Earth, v. 24, 

no. 4, p. 417–444, doi: 10.4294/jpe1952.24.417.
47

Appendix A. Slider-Block Model

In this appendix we explain why we must use a three-dimensional fault model, rather 

than a two-dimensional model.

We wish to construct an earthquake fault model that incorporates concepts from 

thermodynamics. It is natural to consider slider-block models, because these models have 

obvious similarities to the Ising model. (See Rundle and others, 2003; Abaimov, 2009.) In this 

appendix we show that a slider-block model is inappropriate for our application, because it has 

incorrect energy spectrum.

A slider-block model consists of a rectangular lattice of blocks that can slide in the twodimensional plane (see fig. A–1). Each block is connected by springs to its four nearest 

neighbors. In addition, there is a driving force applied to all the blocks. Because of the springs, 

whenever one block slides, the total force on each neighboring block increases, while the total 

force on the sliding block decreases. In other words, tectonic stress is transferred from the sliding 

block to its four neighboring blocks. (See Smalley and others, 1985; Huang and others, 1998.)

An examination of figure A–1 makes plain the similarity between the slider-block model 

and the Ising model. Both models consist of objects arranged in a lattice, where interactions 

occur only between nearest neighbors.

Our objective is to use the Boltzmann probability distribution to choose a random initial 

stress distribution on the fault surface. Because the Boltzmann probability of a state depends on 

the energy of the state, it is essential that each state be assigned a correct energy value. We now

show that the slider-block model does not assign correct energy values and, therefore, that the 

slider-block model is not acceptable for this application.

Let and be coordinates on the plane. Suppose we impose a sinusoidal displacement 

on the blocks:

 

 (79)

where is the displacement of the block located at coordinates , is an 

amplitude, and and are wave numbers in the two coordinate directions. Then the total 

energy stored in the springs scales as 

 

:

 spring 

  here  

  

 (80)

But if you impose a sinusoidal pattern of slip on a fault embedded in a three-dimensional

elastic medium, then the total energy stored in the elastic medium scales as :

 elastic 

 (81)

The difference can be understood as follows. In the case of the slider-block model, the 

springs are stretched or compressed by an amount proportional to , hence the energy in the 

springs is proportional to 

 

. In the case of a fault in a three-dimensional elastic medium, the 

medium immediately adjacent to the fault is distorted by an amount proportional to , hence 

the elastic energy density is proportional to 

 

. But the distortion penetrates into the threedimensional medium, perpendicular to the fault, by a distance proportional to the wavelength, 

which is . So the amount of material subject to the distortion is proportional to , hence 

the total elastic energy is proportional to 

 

 .
48

Because the variation of energy with wavenumber is not correct, we cannot employ a 

slider-block model.
49

Appendix B. Properties of the Basis Displacement Fields

In this appendix we establish the mathematical properties of the basis displacement 

fields, which were defined in equations 8–10.

We use a three-dimensional coordinate system where the direction is 

vertical, the direction is along-strike, and the direction is perpendicular to the fault. The 

fault is the surface . We impose periodic boundary conditions in the and directions 

(but not in the direction), with periods and respectively.

A displacement is represented by a vector . Our basis displacement fields are 

defined as:

 

   

 | |

 ( 

 )  (82)

 

   

 | |

 ( 

 )  (83)

 

 

 

( 

  ) 

 | |

 ( 

 ) (84)

Here 

 and  are arbitrary amplitudes,  is an arbitrary phase, and 

are wavenumbers in the and directions, and is the scalar wavenumber defined as:

 √ 

  (85)

The wavenumbers are assumed to be compatible with the periodic boundary conditions. 

In the main report (equations 8–10), we defined 

 as follows:

 

  

 

 

 

 

 (86)

We will see below that this choice of causes the shear stress to be oriented along-strike. 

For most of this appendix, we allow to be arbitrary.

For the sake of brevity, we introduce the following abbreviations:

 (87)

 (88)

 (89)

 | 

| (90)

 

  (91)

 

  (92)
50

 

    (93)

With these abbreviations, the basis displacement fields can be written as:

 

 

 

 (94)

 

 

 

 (95)

 

 

 

 

 

 (96)

Fault Slip and Exponential Decay

The basis displacement fields are continuous everywhere except at the fault plane , 

where there is a discontinuity due to the factor . The fault slip vector is 

defined as:

   (97)

   (98)

   (99)

We see that the basis displacement field describes a sinusoidal pattern of fault slip, with 

wavenumbers and . The rake angle is determined by the relative size of the amplitudes 

and .

The basis displacement fields depend on only through the factor 

  | |

. So we 

see that the fields decay exponentially with increasing distance, | 

|, from the fault surface, over 

a characteristic distance, .

Derivatives of the Basis Displacement Fields

We need to compute the derivatives of the basis displacement fields. First, note that the 

first derivatives of and are (where we use a subscript comma to denote derivative):

 (100)

 (101)

 (102)

 (103)

 

 (104)
51

Now we can write the first derivatives of the basis displacement fields:

 

 

 

 (105)

 

 

 

 (106)

 

 

 (107)

 

 

 

 (108)

 

 

 

 (109)

 

 

 (110)

 

 

 

 

 

 (111)

 

 

 

 

 

 (112)

 

 

 

 (113)

From the first derivatives, we see that

 (114)

which means that the basis displacement fields are volume-preserving. It follows immediately 

that

 for all (115)

We can now write the second derivatives of the basis displacement fields:

 

 

 

 

 (116)

 

 

 

 

 (117)

 

  

 

 (118)

 

 

 

 

 (119)

 

 

 

 

 (120)

 

  

 

 (121)
52

 

 

 

 

 

 

 (122)

 

 

 

 

 

 

 (123)

 

 

 (124)

From the second derivatives, we see that

 for all (125)

which means that the basis displacement fields satisfy the Laplace equation.

Elastostatic Equation

We show that the basis displacement fields satisfy the elastostatic equation. The strain 

tensor is

 

 

 

( ) (126)

The stress tensor is

 

(127)

where is the shear modulus and is the volumetric Lame parameter. Differentiating, we get

 (128)

We have already shown that each of the terms on the right-hand side is equal to zero. So

 (129)

which shows that the basis displacement fields satisfy the elastostatic equation.

Elastic Energy

Here we calculate the elastic energy associated with the basis displacement field. The 

elastic energy density is:

 

 

 

 

 

 

 

 

 

 

 

 

 

 ( )

 

 

 

 

 

 

 

  

(130)
53

To obtain the elastic energy, we must integrate over all space. For our basis 

displacement fields, the first term on the right-hand side vanishes because we have shown that 

 . To evaluate the second term on the right-hand side, we group expressions as follows:

 

   

   

(131)

 

   

   

(132)

 

   

  

 (133)

Summing the above three equations, we get

 

     

 

 

 (134)

To evaluate the third term in the expression for , we group expressions as follows:

 

 

  

  (135)

 

  

  

  (136)

 

  

  (137)

Summing the above three equations, we get

 

  

 (138)

Combining the above, the elastic energy density is:

 

 

 

 

     

 

 

 (139)

To obtain the total elastic energy, we must integrate over all space. First, note that is 

independent of and . Because we impose periodic boundary conditions in the and 

directions, we have

∫ 

(140)

∫ (141)

Also,

∫ 

 

 

 

 ∫ 

 | | 

 

 

 

 

 

 (142)
54

Therefore the total elastic energy is:

 

 ∫  

 

 

 

 ( (

 

 

 

 

 

 )

 

 

  ) 

(143)

For the special case where shear stress is oriented along-strike, we substitute 

 

 

 

 to get the formula given in the main report (equation 15):

 

 

 

 

 (

 

 

 

 

 

)

 

 

 

 (144)

The above equation is derived using the following algebraic identity:

 ( 

 

 

 

 

 

 

 

 

)

 

 

 

  

 

 

 

 

 

 

 ( 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

)

 

 

 

  

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

  

  

  

 

 

 

 

 

 

 

 

 

 

 

 

 

 

(145)

Orthogonality

We now show that the basis displacement fields are orthogonal. In other words, if we 

form the sum of the basis displacement fields, then the total elastic energy of the sum equals the 

sum of the elastic energies of the individual displacement fields.

Let 

 be the sum of the basis displacement fields:

 

total ∑ 

 

 

(146)

The sum runs over all wave vectors that are compatible with the periodic boundary 

conditions, but with the following two special considerations:

 Because wave vectors and represent the same physical state, the 

sum includes only one member of each such pair of wave vectors.

 Because the zero wave vector 

 represents a state with zero energy 

and zero stress, the sum does not include the zero wave vector.

Using the fact that 

 , we can write the total elastic energy density as:
55

 total 

 

 

 

total total 

 

 

 

total total

 ∑ ∑ 

 

 ( 

   

 

 

 

    

 

 

 

 )

 

 

 

 

 

(147)

The above expression is a pair of sums over , , 

 

, 

 

, , and . Each term in each sum has 

one of the following four functional forms:

 

 

 

 

 

   

 

 

 (148)

 

 

 

 

 

   

 

 

 (149)

 

 

 

 

 

   

 

 

 (150)

 

 

 

 

 

   

 

 

 (151)

If 

 

 

 

 then each of the four expressions above integrates to zero, when integrated 

with respect to and . Therefore, when we integrate over all space, we only need to 

include the terms in the sum that have 

 

 

 

 

 . So the total elastic energy is:

 total ∫ total 

 

 ∑ ∫

 

 

 ( 

        ) 

 

 ∑ ∫ 

 

 

 ∑ 

 

 

(152)

This proves orthogonality.

Normal Stress

We define the normal stress to be the compressive stress perpendicular to the fault 

plane:

 

 

 

 (153)

where we have used the fact that .

We see that the normal stress is not continuous at the fault surface . This shows 

that the basis displacement fields are only an approximate solution to the elastostatic fault slip 

problem.

Note: In the special case, , the normal stress is continuous at the fault 

surface, and so the basis displacement fields are an exact solution to the elastostatic fault slip 
56

problem. The condition holds if the wave vector is perpendicular to the slip 

vector.

Note: We make this approximation because we are not aware of any analytic formulas 

similar to equations 82–84 that give the displacement field throughout space, for the exact 

solution of the elastostatic fault slip problem with a sinusoidal slip pattern.

Note: One can potentially perform a similar calculation to the one given here, without 

making this approximation, despite the lack of analytic solutions. Given a sinusoidal pattern of 

fault slip, the resulting shear stress on the fault plane could be calculated using integral formulas 

similar to those in Andrews (1974) (although it should be noted that the integrals in Andrews 

(1974) are divergent for sinusoidal fault slip). Then, total elastic energy could be computed by 

integrating the dot product of slip and shear stress over the fault plane.

Shear Stress

Let and be the components of shear stress in the direction (along-dip) and 

direction (along-strike), respectively. Their values are:

 ( )

 (

 

 

 

 

 

 

 

 

 

 ) 

 

 

(154)

 ( )

 (

 

 

 

 

 

 

 

 

 

 ) 

 

 

(155)

We see that the shear stress components and are continuous at the fault surface 

 , as stated in the main report.

From here on we are only concerned with the values of shear stress at the fault surface, so 

we drop the exponential term 

 

, which equals when .

If , then the shear stress is oriented along-strike. From the above formulas, we see 

this happens when

 

 

 

 

 

 (156)

Notice that the rake angle (the angle between the horizontal and the shear stress) can be set to 

any desired angle by adjusting the values of and . To prove this fact, note that given any 

numbers and , it is always possible to find values of and that satisfy the simultaneous 

linear equations
57

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

(157)

because of the non-vanishing of the determinant

|

|

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

|

| 

 

 (158)

Let be the two-dimensional shear stress. Then

 

  

 

 

 ( (

 

 

 

 

 

 )

 

 

  )

 

 

(159)

where the sign of can be chosen arbitrarily. To derive the above formula, one needs to use the 

algebraic identity

( 

 

 

 )

 

 

 

 

 

 

 

 

 

     

 

(160)

which can be verified by making the substitution 

 

  

and then expanding both sides of 

the identity.

We previously showed that the elastic energy is

 

 

 

 ( (

 

 

 

 

 

 )

 

 

  ) (161)

Therefore the shear stress is

 √

 

 

 (162)

Notice that the relation between the elastic energy and the shear stress is independent 

of the rake angle. Because the Boltzmann probability distribution depends on energy, this 

implies that the Boltzmann probability distribution cannot provide any information about the 
58

rake angle. So, we must impose a rake angle, such as by requiring the shear stress to be oriented 

along-strike.

For the special case where shear stress is oriented along-strike, we fix the sign of by 

requiring that . Substituting 

 

 

 into the formula for , we get

 ( 

 

  

 

  

 

 

 

 

 

 

 

 

) 

 

 

 

 

 

 

 

(163)

We see that for to have the same sign as , we must have

 √

 

 

 (164)

which is equation 20 in the main report.

The mean-square value of the two-dimensional shear stress , which we denote 

 

, is:

 rms

 

 

 

∫( 

 )

 

 

 

 

(165)

which is equation 21 in the main report.

Total Shear Stress and Completeness

In this part, we confine our attention to the case where shear stress is oriented alongstrike.

As before, let 

 be the sum of the basis displacement fields:

 

total ∑ 

 

 

(166)

The sum runs over all wave vectors that are compatible with the periodic boundary 

conditions, with the special considerations described earlier. The shear stress produced by 

 

is:

 total 

 ∑  

 

 ∑ √

 

 

 

 

 

 

(167)
59

This equation has the form of a two-dimensional Fourier transform. It implies that the 

basis displacement fields are complete in the following sense: Given any pattern of shear stress 

on the fault surface, oriented along-strike and with a finite mean-square, it is possible to find 

some linear combination of the basis displacement fields that produces that pattern of shear 

stress. The linear combination can be found by performing a Fourier transform on the desired 

pattern of shear stress.

The equation also implies that the basis displacement fields are independent in the 

following sense: Any two different linear combinations of the basis displacement fields produce 

different patterns of shear stress on the fault.

We can apply Parseval’s theorem to compute the mean-square value of , which we 

denote by 

 

:

 total

 

 

 

∫( total 

 )

 

 ∑

 

 

 

 

 ∑ rms

 

  

 

 (168)

In other words, the mean-square value of the total shear stress equals the sum of the 

mean-square values of the individual basis displacement fields. This formula (in integral form) 

appears as equation 22 in the main report.

Finally, we note that given a Fourier expansion, the power spectral density is defined to 

be the square of the magnitude of the Fourier coefficients. We see that 

 

 is the power 

spectral density of . In other words, the mean-square shear stresses of the individual basis 

displacement fields are equal to the power spectral density of the total shear stress.
60

Appendix C. Derivation of the Earthquake Rate Formula

In this appendix we show how to derive the earthquake rate formula from the scaling 

relations.

We start by stating three earthquake scaling relations from Hanks (1977; 1979). We give 

them in a slightly more general form than in the main report.

Scaling relation 1 is the Gutenberg-Richter relation (Gutenberg and Richter, 1954). Let 

be an earthquake magnitude, and let be the rate of occurrence of earthquakes with magnitude 

 . Then

 (169)

where and are some constants.

Scaling relation 2 is the relation between magnitude and seismic moment (Hanks and 

Kanamori, 1979). Let be seismic moment, which is defined as

 

 

(170)

where is fault slip, is the radius of the fault rupture area, and is some constant. It is 

conventional to take where is the shear modulus, but this particular choice of is not 

necessary for the present discussion. Then

 (171)

where and are some constants.

Scaling relation 3 is the statement that earthquakes are self-similar or scale-invariant. It is

 (172)

where is some constant.

Substituting equation 172 into equation 170, we have

 

 

 (173)

Substituting this into equation 171 and solving for , we get

 

 

 

 

 

 

 

 

 (174)

Substituting this into equation 169, we get

 ( 

 

 

) 

 

 

 

 

 (175)

Differentiate the above formula with respect to to obtain
61

 

 

 

 

 

 

 (176)

The above differential equation for as a function of has solution

 

 

(177)

 

 

 

 

 

 (178)

where is some constant. Note that because the earthquake rate, , decreases as the 

radius increases.

Substituting in the empirical values and , we obtain equation 42 in the main

report:

 

 

 

 

 

 (179)
62

Appendix D. Derivation of the Power Spectral Density

In the main report (equation 47) we demonstrated a scaling relation for stress 

perturbations on the fault surface. We showed that 

 (180)

where denotes the total area of the fault surface occupied by positive stress 

perturbations of radius .

In this appendix we derive the power spectral density of stress, which is consistent with 

the above scaling relation. Our end result will be a formula that gives the asymptotic behavior of 

the power spectral density, in the limit of large wavenumber.

Notations and Definitions

First we recall some notations from the main report. The fault plane has coordinates 

and , where runs vertically and runs along-strike. We impose periodic boundary 

conditions, with periods and , so that the fault area is . We use and to denote 

wavenumbers in the two coordinate directions, and the scalar wavenumber is √ 

  

.

We denote the shear stress pattern on the fault as 

 . We think of it as being the 

sum of plane waves (essentially the same thing as a Fourier transform):

 total 

 

 

 

 ∫ 

 (181)

The constant 

 

is called the density of states and is explained in the main report 

(equation 18); its value is not important for the present discussion. The individual plane waves 

take the form

 

 (  ) (182)

where and 

 are random amplitudes and phases, respectively.

Let 

 

  be the mean-square value of the individual plane wave. In the above 

notation, it is ( 

 )

 

 . Then 

 

  is the (two-dimensional) power spectral density 

of the stress pattern.

The shear stress distribution on the fault is generated by some random process. Whenever 

 e use the ord “expected” or enclose an expression in angle brackets, e mean the average 

value obtained over many such randomly generated stress distributions. For example, 

 

 

  denotes the power spectral density for some particular stress distribution, while 

〈 

 

  〉 denotes its expected or average value over many stress distributions.

We assume that the expected power spectral density 〈 

 

  〉 is isotropic, that is, it 

depends on the scalar wavenumber but not on the direction of the wave vector. Note that the

power spectral density 

 

  for a single stress distribution is not isotropic. Only its 

expected value 〈 

 

  〉 is isotropic. This, in part, is why we must work with expected 

values and not individual stress distributions.

The expected mean-square value of the shear stress is defined to be
63

⟨ total

 

⟩ ⟨

 

 

∫( total 

 )

 

 ⟩ (183)

Using Parseval’s theorem, this can be ritten as

⟨ total

 

⟩ ⟨

 

 

 ∫ rms

 

 ⟩ 

 

 

∫ 〈 rms

 

  〉 

 

 

(184)

where the final integral depends on the fact that 〈 

 

  〉 is isotropic. The constants 

appearing in front of the integrals are explained in the main report (equations 18 and 19); their 

values are not important for the present discussion.

Filtered Stress Fields

Next we introduce the concept of a filtered stress field, which contains only the highfrequency components of the stress distribution. This is useful because the high-frequency 

components of the stress distribution are responsible for the small-scale stress perturbations, 

which enter into the scaling relation. Also, we are trying to find the asymptotic behavior of the 

expected power spectral density, in the limit of large wavenumber.

Let denote the portion of the stress field that is due to wavenumbers 

 . In other words, is the result of applying a high-pass spatial filter to the 

stress distribution 

 , chosen so that the filter only passes wavenumbers . If 

 

 is constructed by forming the sum of plane waves with various wave vectors, then 

 is constructed by forming the sum of the same plane waves, except that only 

those plane waves with wavenumbers are included in the sum. Using the notations 

introduced above,

 total 

 

 

 ∫ 

 

 

  

 (185)

The expected mean-square value of the filtered shear stress is:

⟨ total

 

 ⟩ 

 

 

∫ 〈 rms

 

  〉 

 

 

 (186)

Notice that the limits of integration run from to because the filtered stress only includes 

wavenumbers .

Expected One-Point Statistics of Filtered Stress

We need to know more about the filtered stress fields than the expected mean-square 

value. To do this, we study their expected one-point statistics, which is the expected statistical 

distribution of stress values at individual points on the fault surface.
64

For each , define to be a real-valued random variable which has the expected 

one-point statistics of the filtered stress field . What this definition means is 

that the cumulative density function of is defined to be:

 ⟨

 

 

 total ⟩ (187)

where is the probability that a randomly selected value of is greater than or equal 

to , and where denotes the area of the portion of the fault surface 

where the filtered stress field is greater than or equal to .

In other words, the probability that is greater than or equal to is equal to the 

expected fraction of the fault surface where the filtered stress is greater than 

or equal to .

With this definition, the mean of is zero, and the variance of is the expected meansquare stress:

 (188)

 ⟨ total

 

 ⟩ 

 

 

∫ 〈 rms

 

  〉 

 

 

 (189)

The proof of these formulas is lengthy. A subsection at the end of this appendix shows 

how to prove these formulas, and also addresses some technical issues related to the definition of 

 .

Small-Scale Stress Perturbations

Now we need to make the connection between and the total area of small-scale 

positive stress perturbations, which is what appears in our scaling relation.

Because contains all wavelengths , we say that 

 contains all stress perturbations of radius for some constant . (It 

would be reasonable to take , but, as we will see, the value of has no effect on our final 

result.)

We say that a positive stress perturbation occurs when the stress 

equals or exceeds some constant level . Then, the total area of all positive stress perturbations 

of radius is

 total (190)

Notice that the threshold stress level is independent of . This is a key assumption. It 

comes from assumption of earthquake self-similarity, which states that the stress drop in an 

earthquake is independent of the size of the rupture area.

The expected total area of all positive stress perturbations of radius is

〈 〉 〈 total 〉 (191)
65

As varies, we require the expected area to vary in accordance with our scaling relation

〈 〉 (192)

Combining the above two equations, we obtain the scaling relation

 (193)

where the constants , , and have been absorbed into the implied constant of 

proportionality. This is the scaling relation we use in the calculations to follow.

Method of Calculation

Here we explain the method we will use to start with the scaling relation above, and 

calculate the expected power spectral density.

For convenience, let us define by

 (194)

Then, recalling that the expected power spectral density 〈 

 

  〉 is isotropic, we have

 

 

 

∫ 〈 rms

 

  〉 

 

 

(195)

 

 

 

 

〈 rms

 

  〉 (196)

and so the expected power spectral density is:

〈 rms

 

  〉 

 

 

 

 (197)

Our objective is to find asymptotically in the limit . Then, the above 

formula gives an asymptotic expression for the expected power spectral density 〈 

 

  〉.

Let us rewrite our scaling relation with an explicit constant of proportionality :

 (198)

Notice that has dimensions of inverse length, because the wavenumber has dimensions of 

inverse length.

The remainder of the calculation proceeds as follows. First, we make an assumption 

about the probability density function of the random variables . Then, we start from the 

scaling relation and derive a formula for or . If it is not possible to derive an exact 

formula in closed form, then we derive a formula that is asymptotically true in the limit . 
66

Then, we immediately obtain an asymptotic formula for the expected power spectral density 

〈 

 

  〉.

We will first assume that each is a Gaussian random variable. This assumption is the 

most physically reasonable, and it yields the formula given in the main report (equation 48). 

Then, to explore the sensitivity of our results to the assumed probability distribution of , we 

repeat the calculation using three other probability distributions: the Laplace, exponential, and 

Pareto distributions. We will see that the Gaussian, Laplace, and exponential distributions, which 

have tails that fall off exponentially, all produce very similar results: the power spectral density 

falls off as 

 modified by a logarithmic term. For the Pareto distribution, which has a tail that 

falls off as a power law, the power spectral density falls off as a power law that is faster than 

 

 

, but slower than 

 

. The implication of these results, especially the three with 

exponential tails, is that the power spectral density is not very sensitive to the precise probability 

distribution of the random variable .

Before beginning, a word of warning about the manipulation of asymptotic expressions is 

needed. We say that functions and are asymptotic in the limit , and we write 

 , if and only if, 

  in the limit . Some operations that are valid 

on equalities are not valid on asymptotic relations. For example, if , then it follows 

that , but not that . In particular, it is not generally

valid to differentiate an asymptotic relation; from , it does not follow that 

 

 

 

 

 .

Calculation Assuming Gaussian Distribution

Here we assume that each is a Gaussian random variable, with zero mean and 

variance .

We begin with some preliminaries. The complementary error function is defined as

 

 

√ 

∫ 

 

 

 

 

 

 (199)

Its derivative is

 

 

 

√ 

 

 

 

 (200)

The complementary error function, and its inverse, have asymptotic formulas (National Institute 

of Standards and Technology, 2011):

 

 

√ 

 

 

 

as (201)

 ( 

 

 

 

 

 )

 

as 

(202)

For our purposes, we use a simpler (albeit more slowly converging) asymptotic formula 

for the inverse complementary error function,
67

 ( (

 

√ 

))

 

as (203)

This is also a valid asymptotic formula because

 

 

 

 

 

 

 (204)

This completes the preliminaries. If is Gaussian with zero mean and variance , then 

our scaling relation

 (205)

Becomes

 

 

 

 

 

 (

 

√ 

) (206)

Differentiate the above equation with respect to , using the formula for the derivative of 

the complementary error function, to obtain

 

 

 

 

 

 

(

 

√ 

)

 

(

√ 

√ 

 

 

 )

 

 (207)

Noting that as , use the asymptotic formula for the complementary error 

function to get

 

 

 (

 

√ 

) 

√ 

√ 

 

 

 

as (208)

Using the asymptotic formula for the inverse complementary error function,

 

√ 

 (

 

 

) 

 

as (209)

where √ is a parameter with dimensions of length.

(From a strictly mathematical standpoint, one could eliminate and write in place 

of , and still have a correct asymptotic formula. We retain because the wavenumber 

has dimensions of inverse length, while the logarithm function requires a dimensionless 

argument, so it would be physically incorrect to omit .)

Combining the three previous formulas,

 

 

 

 

 

 

 

 

 (

 

 

)

 

 as (210)
68

which can be re-written as

 

 

 

 

 

 

 

 

 

as (211)

So the expected power spectral density is

〈 rms

 

  〉 

 

 

 

 

 

 

 

 

as (212)

This is equation 48 in the main report.

Remark: It is interesting to note that, if it were permitted to differentiate an asymptotic 

formula, the above result could be obtained much more quickly. Using the asymptotic formula 

for the inverse complementary error function,

 

 

 

 

( (

 

 

))

 

 

 

 

 

 

 

as 

(213)

Differentiating the right-hand side of the above formula with respect to , were it permissible to 

do so, would immediately yield the same asymptotic result for .

Calculation Assuming Laplace Distribution

Here we assume that each is a Laplace random variable, with zero mean and variance 

 . The Laplace distribution is a symmetric double-tailed distribution, with probability 

density proportional to | | . By comparison, the Gaussian distribution is a symmetric 

double-tailed distribution with a probability density proportional to 

 

 

 

 . So, this 

calculation explores the effect of using a distribution that falls off more slowly than Gaussian.

The Laplace distribution is defined to have the cumulative density function

 

 

 

 

 √ 

assuming (214)

The restriction on is not a problem, because is the threshold stress for defining positive stress 

perturbations, so physically we will always have .

Our scaling relation becomes

 

 

 

 

 

 

 √ 

 (215)

The above formula can be solved for , obtaining

 

 

 

 

 

(216)
69

where is a parameter with dimensions of length, and where we assume . 

Differentiating the above equation, we get

 

 

 

 

 

 

 

 

 

 (217)

So the expected power spectral density is

〈 rms

 

  〉 

 

 

 

 

 

 

 

 

 (218)

Notice that this agrees with the Gaussian formula, except that it has 

 

in place of 

 

 

. The fact that the two formulas are so similar helps give confidence in our use of the 

Gaussian formula.

Calculation Assuming Exponential Distribution

Here we assume that each is an exponential random variable, with zero mean and 

variance . The exponential distribution is a single-tailed distribution, with probability 

density proportional to . By comparison, the Gaussian distribution is a symmetric 

double-tailed distribution with a probability density proportional to 

 

 

 

 . So this 

calculation explores the effect of using a distribution that falls off more slowly than Gaussian, as 

well as the effect of using an asymmetric distribution. (Physically, an asymmetric distribution 

might be interpreted as saying that the stress cannot be much lower than its mean value, but can 

be much higher than its mean value.)

The exponential distribution (shifted to have zero mean) is defined to have the 

cumulative density function

 

 √ 

assuming √ (219)

The restriction on is not a problem, because is the threshold stress for defining positive stress 

perturbations, so physically we will always have .

Our scaling relation becomes

 

 

 

 √ 

 (220)

The above formula can be solved for , obtaining

 

 

 

 

 

(221)

where is a parameter with dimensions of length, and where we assume . 

Differentiating the above equation, we get
70

 

 

 

 

 

 

 

 

 

 (222)

So the expected power spectral density is

〈 rms

 

  〉 

 

 

 

 

 

 

 

 

 (223)

This is the same formula obtained with the Laplace distribution.

Calculation Assuming Pareto Distribution

Here we assume that each is a Pareto random variable, with zero mean and variance 

 . The Pareto distribution is a single-tailed distribution, with probability density 

proportional to a power law 

 

. We require for the Pareto distribution to have finite 

variance. So, this calculation explores the effect of using a distribution that falls off much more 

slowly than Gaussian.

The Pareto distribution (shifted to have zero mean) is defined to have the cumulative 

density function

 [√

 

 

 

√ 

√ √ 

]

 

 ith √

 

 

 

√ 

(224)

The restriction on is not a problem, because is the threshold stress for defining positive stress 

perturbations, so physically we will always have .

Our scaling relation becomes

 

 

 [√

 

 

 

√ 

√ √ 

]

 

 (225)

The above formula can be solved for , obtaining

 

 

[(

 

 

 

)

 

 ]

 

(226)

where we define

 

 

 

√ 

(227)
71

 

 (

 

 

)

 

 (228)

Differentiating the above equation, we get

 

 

 

 

 

 

 

 

[(

 

 

)

 

 ]

 

(

 

 

)

 

(229)

and so

 

 

 

 

 

 

 

 

(

 

 

)

 

as (230)

So the expected power spectral density is

〈 rms

 

  〉 

 

 

 

 

 

 

 

as (231)

Recalling that we need for the Pareto distribution to have finite variance, we see 

that the exponent satisfies . So the asymptotic falloff of the expected power 

spectral density is a power law that is faster than 

 

, but slower than 

 

.

Technical Issues for One-Point Statistics

Here we discuss some technical issues relating to the definition of and the formulas 

for its mean and variance. Our goals are to (1) define the concept of a random variable that has 

the expected one-point statistics of the filtered stress distribution; (2) show that its mean is zero; 

and (3) show that its variance equals the expected mean-square filtered stress.

Removing Infinities

We introduced periodic boundary conditions to make the fault area finite. This avoids a 

number of mathematical and conceptual difficulties. For example, a finite fault makes it 

straightforward to talk about the mean-square value of the stress or the total elastic energy of a 

displacement field.

Another advantage of periodic boundary conditions is that it makes the wave vector 

assume discrete values. But there are still an infinite number of wave vectors. To pick a stress 

distribution at random, one must (in theory) pick an infinite number of wave amplitudes. This 

creates some difficulties. For example, the resulting stress distribution would not be a continuous 

function (it would instead be a so-called 

 

function); it is not straightforward to define the value 

of the stress at some given point on the fault; and it is difficult to work with an infinite number of 

random variables.

We can remove the infinity by imposing a maximum wavenumber . That is, we 

construct our stress field by summing only plane waves with | 

| and | | . The 

limit is treated as arbitrarily large, much larger than any wavenumber of interest. With this 
72

constraint, there are only a finite number of admissible wave vectors , so a stress 

distribution can be chosen by selecting a finite number of wave amplitudes. As a result, the set of 

possible stress fields is a finite-dimensional vector space, the stress fields are continuous, and the 

value of stress at a given point is completely well-defined.

The maximum wavenumber is analogous to the fault dimensions and , in that 

it is a large number introduced to eliminate infinities from the calculation, with the expectation 

that the final result does not depend on its value so long as the value is very large.

(Technical remark: We impose the conditions | 

| and | | , instead of 

the condition , because the former condition guarantees that two-dimensional Fourier 

transforms converge pointwise in the limit (Fefferman, 1971), whereas it is unknown 

if the latter condition does so. For our purposes, the exact form of the wavenumber cutoff is not 

important.)

Stress Distributions and Wave Amplitudes

Let 

 denote the shear stress on the fault, which we presume is generated by 

some random process. Due to the existence of a maximum wavenumber, the stress can be 

expressed as a finite sum of cosines and sines:

 total 

 ∑  

 

 ∑ 

 

 

 

(232)

Recall from the main report (equation 18) that we do not include any terms with zero 

wave vector, 

 , because a uniform displacement field generates no stress and no 

elastic energy. Physically, this means that 

 is the deviation in stress from some 

unspecified mean stress value.

Let denote the set of possible wave amplitudes . Then is a vector space of 

finite dimension , and each point in corresponds to a stress distribution according to the 

above formula.

To randomly select a stress distribution, all that is needed is to choose random numbers 

 , that is, to choose a point in . We assume that the random selection is made in 

accordance with a probability density function , which is a non-negative measurable 

function with an integral of .

(Informally, a function is called measurable if it is sufficiently well-behaved so that it can 

be integrated. Rigorous mathematical definitions may be found in, for example, Royden (1968).)

If is a real-valued measurable function on and, hence, also a real-valued function on 

the set of possible stress distributions, then the expected value of is

〈 〉 ∫  

 

 (233)

The subscript emphasizes that the integral is over the set of wave amplitudes, or stress fields.

Recall that denotes the portion of the stress field due to 

wavenumbers . In other words, it is the result of applying a high-pass filter to 

 , 

where the filter only passes wavenumbers . (Of course, we assume .) Then 
73

 has the same expansion as a sum of cosines and sines, with the same 

amplitudes, except omitting any terms with 

  

 

. So, we have

 total ∑ 

   

 

  

 ∑ 

 

 

  

 

(234)

Fault Surface and Mean Values

Let denote the fault surface. Then consists of points , subject to the periodic 

boundary conditions. The fault surface has a finite area, which is .

By Parseval’s theorem, the mean-square value of the stress field 

 equals onehalf the sum of the wave amplitudes:

 total

 

 

 

∫( total

  )

 

 

 

 

 

 

∑( 

  )

 

 

 

 

 

∑ ( 

  )

 

 

 

(235)

The subscript emphasizes that the integral is taken over the fault surface . Likewise, the 

mean-square value of the filtered stress field is

 total

 

 

 

 

∫( total )

 

 

 

 

 

 

∑ ( 

  )

 

 

  

 

 

 

∑ ( 

  )

 

 

  

 

(236)

We see that applying the high-pass filter reduces (or leaves unchanged) the mean-square value:

 total

 

 total

 

 (237)

If is a real number, define the step function to be

 

 {

 if 

 if 

 (238)

Then is a measurable function (although it is not continuous). If is a measurable 

real-valued function on , then is also a measurable function. We define:

 

 ∫  

 

 (239)
74

We say that 

 is the area of the portion of the fault surface where the condition 

 

 is true. That is, it is the area of the portion of the fault surface where the function 

value is greater than or equal to .

Product Space

Let denote the product space of and . By definition, an element of is a 

tuple where and . Then is a space with 

dimensions.

The function 

 is a probability density function on , as shown by:

 

 

∫ 

 

 

 

 

 

∫ [∫ 

 

] 

 

 

 ∫ 

 

 

 

(240)

The subscript emphasizes that the first integral is taken over the product space . It 

can be evaluated as an iterated integral by Tonelli’s Theorem ( oyden, theorem 12.20). electing 

a random point according to the probability density function 

 

 is equivalent to selecting a random point according to the uniform probability 

distribution, and independently selecting a random set of wave amplitudes 

according to the probability density function .

Definition of Expected One-Point Statistics

Here we define a random variable , which has the expected one-point statistics of 

 .

Conceptually, choosing a value of the random variable could be done like this: (1) 

Select a random point on the fault surface, according to a uniform probability 

distribution. (2) Select a random set of wave amplitudes according to the probability 

density function . This is equivalent to selecting a random stress field. (3) Let be the 

corresponding value of , which can be computed from the expansion of 

 into a finite sum of cosines and sines.

We want a mathematically precise definition of . For , let be a real-valued 

function on the product space , defined as follows:

 

 total (241)

Then is a continuous function and, hence, also a measurable function, because 

 is the sum of a finite number of cosine waves and sine waves, all of which 

are continuous functions.

We define the random variable to have the probability measure that is induced by the 

function . That is, we define

 

 

 

∫ 

 

 

(242)
75

where is the probability that is greater than the real number . 

By the Lebesgue Convergence Theorem (Royden, theorem 11.16) we have

 

 

 (243)

 

 

 (244)

which shows that the definition does indeed describe a probability distribution. By Fubini’s 

Theorem (Royden, theorem 12.19), the integral over can be evaluated using iterated 

integration:

 

 

 

∫ [∫ 

 

 

] 

 

 

 

 

 

∫ total 

 

 

 ⟨

 

 

 total ⟩

(245)

which is just the expected fraction of the fault surface where the value of the filtered stress 

 equals or exceeds .

The definition of has the following consequence (Royden, proposition 15.1). Suppose 

that is a measurable function, and assume that

 

 

∫ | 

 | 

 

 (246)

(note the absolute value in the integrand). Then the expected value of with respect to the 

random variable is given by

⟨ ⟩ 

 

 

∫ 

 

 

 (247)

The subscript on the angle bracket denotes the expected value when is chosen at 

random according to the probability distribution of . It is not to be confused with unadorned 

angle brackets, which always denote the expected value when a stress field is chosen at random 

according to the probability density function . Furthermore, by Fubini’s Theorem 

(Royden, theorem 12.19), the integral over can be evaluated using iterated integration:
76

⟨ ⟩ 

 

 

∫ [∫ 

 

 

] 

 

 

 ⟨

 

 

∫ 

 

 

⟩ 

(248)

Mean and Variance of Expected One-Point Statistics

We now compute the mean and variance of .

As discussed in the main report, we assume that the expected mean-square stress is finite. 

This is because the fault has finite strength, and also because infinitely large stresses are not 

physically realistic. So

 ⟨ total

 

⟩ ⟨ total

 

 ⟩ ⟨

 

 

∫( total )

 

 ⟩

 ⟨

 

 

∫( 

 )

 

 

 

⟩

 

 

 

∫ [∫( 

 )

 

 

 

] 

 

 

 

 

 

∫ ( 

 )

 

 

 

 

(249)

In the above equation, it is valid to change the iterated integral into an integral over by 

Tonelli’s Theorem ( oyden, theorem 12.20). We therefore have

⟨ 

 

⟩ ⟨

 

 

∫( 

 )

 

 

 

⟩ ⟨ total

 

 ⟩ (250)

By the Holder Inequality (Royden, theorem 11.25 and problem 11.21), we have

 

 

∫ | 

 |  

 

 [

 

 

∫ ( 

 )

 

 

 

 

]

 

 

(251)

Therefore,
77

⟨ ⟩ ⟨

 

 

∫ 

 

 

⟩

 ⟨

 

 

∫ total 

 

⟩ 

(252)

In the above equation, the final equality to zero is because each by 

construction has zero mean (that is, it consists of cosine and sine waves with no DC component), 

assuming .

Combining the above equations,

 ⟨ ⟩ (253)

 ⟨ 

 

⟩ ⟨ total

 

 ⟩ (254)

These are the results shown earlier in this appendix.
78

Appendix E. Algorithms for Energy Rolloff

In this appendix we present algorithms for working with the energy rolloff function, and 

for selecting random energy values in accordance with the Boltzmann distribution that have the 

desired rolloff.

Energy Rolloff Function

Recall from the main report (equations 31 and 51) that we defined an energy rolloff 

function where is the scalar wavenumber. For each wave vector , there is a basis 

displacement field, and its expected energy is

〈  〉 

 

 

 (255)

We chose the following functional form for :

 

 

 

  

 (256)

The constants and are adjustable parameters with dimensions of length. Generally, we 

take and adjust to obtain the desired correlation length. Recall that the correlation length 

is defined to be ̂, where ̂ is chosen to maximize the value of ( ̂) ̂ . The definition 

implies that ̂ is the solution to the equation 

 

 . Evaluating the derivative, we 

obtain the equation

 ( ̂)

 

[

 

 

( ( ̂))

 

 ( ( ̂))

 ̂

 ̂

] (257)

Given a desired value of the correlation length , one can substitute ̂ and then 

solve the above equation numerically to obtain suitable values of and . The right-hand side of 

the above equation is a monotone increasing function of both and , and so the bisection 

method can be used to find numerical solutions. For example, our computer code lets the user 

impose a constraint, , where and are user-specified values. Substituting 

 into the above equation yields an equation for , which can be solved by the bisection method.

Note that if then we have .

In addition to describing the decrease in energy with increasing wavenumber, the energy 

rolloff function also establishes a relation between the thermodynamic parameter and the 

expected mean-square stress (equation 34). Given one, the other can be calculated using

⟨ total

 

⟩ 

 

 ∫ 

 

 

 

 (258)

Given our formula for , the improper integral can be evaluated numerically by 

splitting it into the sum of two integrals:
79

∫ 

 

 

 

 ∫

 

 

 

  

 

 

 ∫

 

 

 

 

 

    

 

 

 

(259)

The value of should be chosen so that lies in the tail of the function 

 

. In practice, 

we use . The second integral arises from the change of variable 

 

 

. Each of the two integrals on the right-hand side can be evaluated numerically 

using adaptive impson’s rule. Note that the integrand in the second integral is finite at 

despite the appearance of 

 

in the numerator.

State Restriction

For each wave vector , we impose a maximum allowed amplitude 

for the basis displacement field:

 

 max  (260)

Subject to this restriction, the amplitude, , and phase, , are chosen randomly 

according to the Boltzmann probability distribution for a system with two degrees of freedom. 

Each is chosen in such a way so that the expected energy satisfies 〈 〉 

 .

We now develop a method for generating the random amplitudes and phases. Recall the 

relation (equation 15) between energy and displacement amplitude:

 

 

 

 

 

  

 (261)

where we define

 

 (

 

 

 

 

 

)

 

 (262)

We see that imposing a maximum amplitude is equivalent to imposing a 

maximum energy 

   

 

, with the restriction that

 

 max  (263)

For convenience we temporarily write in place of , and similarly for other 

variables. Define

 

 max (264)
80

or equivalently

 max 

 

 

 (265)

Then the partition function is

 ∫ ∫ 

 

 max

 

 

 

 ∫ 

 

 max

 

 

 

 

∫ 

 

 max

 

 

 

 

(266)

We changed variables from to using . Also note that if we define 

 and , then , which justifies the first integral 

above. The cumulative probability distribution for is

 

 

 

 

 

∫ 

 

 

 

 

 

 

 

 

(267)

where denotes the probability that the energy is less than or equal to some value . 

The implication of the above formula is that if we define

 

 

 

 

 (268)

or equivalently

 

 

 

 (269)

then is uniformly distributed between and . This provides a practical method for generating 

a random value of the energy in accordance with the Boltzmann distribution: generate a 

random value uniformly distributed in and apply the above formula to obtain . A 

random value of the phase can be generated as where is uniformly distributed in 

 . The expected energy is obtained by differentiating the partition function (equation 4):
81

〈 〉 

 

 

 

 

 

 

[ 

 

 

 ] (270)

Notice that can range from to , and in the process 〈 〉 ranges from to . So, 

given any desired value of the expected energy 〈 〉 in the range from to , the above 

equation can be solved numerically to obtain the required value of . Recalling that 〈 〉 

 , the equation that must be solved to obtain is:

 

 

 

 (271)

The equation can be solved using Ne ton’s method, except that the starting value must 

be carefully chosen. Because 

 

 as , Ne ton’s method can converge very slowly 

or fail to converge if the starting value is poor. The following technique can be used to select a 

proper starting value. Define . (The value is chosen so that 

 

 , noting 

that 

 

 ). If , then the starting value is:

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 (272)

The above formula contains the first six terms of the power-series expansion of as a function 

of .

Otherwise, if , define an increasing sequence of values by 

  

 

. The starting value for Ne ton’s method is , where is the smallest index such that 

 . If becomes very close to before the condition is satisfied (which 

may occur if is very close to ), then take .

Remark: The starting value works for the following reasons. If , the power 

series gives a fairly accurate approximation of . For the case , note that 

 

 

and 

 for all , hich implies that convergence of Ne ton’s method is guaranteed 

provided that the starting value is greater than the true value of . Convergence can be very slow 

due to the fact that 

 

 as . But Ne ton’s method converges at least as fast as 

bisection if the starting value and the true value of lie within an interval where the derivative 

changes by a factor of (or less). The sequence is constructed so that 

 

 

 

 , 

and the starting value is the first value in the sequence that is greater than or equal to the true .

Remark: As approaches , the recurrence formula 

 

 

suffers a loss 

of precision. This can be avoided by writing , because then the recurrence becomes 

 

 which does not suffer a loss of precision. Note that if we also write then 

the equation to solve is

 

 

 

 (273)

To avoid loss of precision, Ne ton’s method can be used to solve this equation for , and then 

is obtained as .

Remark: In our computer implementation of the method, we impose a maximum value on 

 :
82

 

 

 (274)

If the computed violates the above limit, then is set equal to its maximum allowed value. 

This has the effect of limiting the random energy to a maximum of times its expected 

value. Because the energy follows an exponential distribution, this is equivalent to standard 

deviations above its mean. Note that this limit only comes into play for , and so it 

affects only long-wavelength (small wavenumber) Fourier modes where the energy rolloff 

is close to . If this was not done, then the random stress pattern could be dominated by a small 

number of long-wavelength Fourier modes with extremely large amplitudes. We consider this to 

be a form of discretization error and adopt the above limit to prevent it. (There is nothing special

about the value ; there is a range of values that would work as well.)
83

Appendix F. Mathematics of the Selection Algorithm

In this appendix we explain the mathematics behind the selection algorithm. As explained 

in the main report, the purpose of the selection algorithm is to select the fault surface from a 

much larger distribution of randomly generated stresses.

Recall (from equations 58 and 59) that the large stress distribution is defined on an 

 grid of points :

 

 for (275)

 

 for 

(276)

where and are integers, and where and are the physical size of the grid in the 

and directions. The grid is considered to satisfy periodic boundary conditions. Because the 

points are discrete, a function is nothing more than function values, and 

so may be regarded as an element of an -dimensional vector space. Given two functions,

 and , we can calculate their dot product within this vector space:

 ∑ 

  

 

 

(277)

The sum runs over the possible values of and is nothing more than the usual dot 

product in a slightly unusual notation.

The dot product can be used as a measure of the extent to which two vectors point in the 

same direction. It is an imperfect measure, but in general, a large positive dot product indicates a 

likelihood that the two vectors point in roughly the same direction. Expressed in the language of 

functions, a large positive dot product indicates a likelihood that the two functions are similar.

The Convolution

Recall from the main report (equation 72) that the complex-valued randomly generated 

stress field is

 total

 

  ∑ √

 

 

 

 ( 

 )

 

(278)

where and are wavenumbers in the and directions, are randomly selected 

energies, and are randomly selected phases.

We are given a real-valued selection template that describes the desired pattern 

of stress on the fault surface, under the assumption that the upper left corner of the fault surface 

is aligned at the origin 

 . Then, we form the convolution (equation 73):

 total

 

  ∑  total

 

 

 

 

(279)
84

For each , the sum on the right-hand side is a dot product of two functions. One function 

is the stress, 

 

. The other function is the selection template, , shifted so that the upper left 

corner of the fault surface lies at the point . To see that the selection template is shifted in 

this way, note that for 

 the term appearing in the sum is 

 

 , 

which shows that the origin of the selection template is aligned with point of the stress 

field.

So if we take ̃ ̃ to be the point which maximizes the absolute value of the 

convolution, then ̃ ̃ is the point where we should place the upper left corner of the fault 

surface so as to maximize the absolute value of the dot product of the stress field with the 

selection template.

The Phase Shift

Recall from the main article (equation 74) that is the complex phase of the convolution,

 

 

 ̃ ̃  , at the point where its absolute value is maximized, that is,

 total

 

 ̃ ̃  | total

 

 ̃ ̃  | 

 

 (280)

Now define

 

 

  total

 

  

 

 (281)

Then

 

 

 ̃ ̃  total

 

 ̃ ̃  

 | total

 

 ̃ ̃  | (282)

In other words, 

 

 ̃ ̃  is real and positive. In fact, it is the largest positive value 

that can be obtained by considering any point in the stress field and any phase. Now define

 

 ( 

 

  )

 ∑ √

 

 

 

 

 

 

(283)

which is the phase-shift formula from the main report (equation 75). Because the selection 

template is real, we have

 ̃ ̃  ( 

 

 ̃ ̃  ) | total

 

 ̃ ̃  | (284)

So we have also maximized the value of . This means that the phase-shifted stress 

field 

 is likely to be similar to the location-shifted selection template, with the shift 

chosen to place the upper left corner of the fault surface at the point ̃ ̃ . Indeed, the chosen 

point ̃ ̃ and phase are likely to maximize the degree of similarity.
85

Appendix G. Post-Processing

In this appendix we describe post-processing steps that can be used to prepare the random 

stress distribution for use in a dynamic rupture simulation.

The main report describes how to construct a raw random stress field on a fault surface. 

Before the stress field can be used as input to a dynamic rupture simulation, it is generally 

necessary to perform some post-processing steps. The steps required will depend on the details 

of the model, such as the friction parameters, the variation of normal stress with depth, and the 

characteristics of the simulation software. Here we describe typical post-processing steps that 

may be used, in the order that they are applied.

Stress Field Rotation

Recall that the random stress distribution on the fault surface is cut out from a much 

larger stress distribution. The cut-out portion is selected by finding the portion of the large stress 

distribution that is most similar to a user-specified selection template , where and 

are the along-dip and along-strike coordinates on the fault surface. See the main report (equation 

73) for details.

If the selection template has symmetry under rotation by 180 degrees, as is the 

case for the selection template shown in figure 8, then it is possible to rotate the stress field by 

180 degrees around the center of symmetry of the selection template. The effect of performing 

such a rotation is that the stress pattern, which is originally near the Earth’s surface, is moved to 

the bottom of the fault and vice versa. (Note that if the center of symmetry of the selection 

template does not coincide with the center of symmetry of the fault surface, then performing the 

rotation also causes a slightly different portion of the large stress distribution to be used when the 

fault surface is cut out.)

For many dynamic rupture simulations, particularly in cases where the normal stress 

increases with depth, high stress near the bottom of the fault is more effective at driving the 

rupture than high stress near the Earth’s surface. o, one can perform a test to check if the 

original stress field has more high stress in the upper part of the fault than in the lower part and, 

if so, rotate the stress field by 180 degrees.

To test whether or not to perform the rotation, we calculate

∑ 

 

 

(285)

where the sum runs over all points on the fault surface, is a user-specified 

function called the rotation template, and is the random shear stress. The above sum is 

calculated for both the rotated and un-rotated stress fields, and whichever one has a larger sum is 

used. The rotation template is typically chosen to be zero at the Earth’s surface, and 

positive near the bottom of the fault, so that stresses at depth contribute more to the sum than 

stresses near the Earth’s surface.

Filtering

Simulation codes typically cannot resolve details smaller than 3 to 5 grid points. Our 

randomly generated stress field contains details down to the level of individual grid points. It is 
86

desirable to suppress such details, so as not to exceed the resolution of the simulation code. We 

do this by applying a weighted moving-average filter to the stress field:

 

 ∑ ∑   

 

 

 

 

(286)

where and are the spacing between adjacent points on the fault surface. The coefficients 

 are called the filter kernel. We take them to be bilinear, which is defined as:

 

 | | | | (287)

The constant is chosen so that

∑ ∑ 

 

 

 

 

 

 (288)

But the choice of filter kernel is not important. The moving-average filter acts as a lowpass spatial filter. Although it can be evaluated using Fourier transforms, it is typically better to 

just evaluate the summation directly, because typically and are quite small. For the 

examples included with this report, we used .

An issue with a moving-average filter is what to do near the edges of the fault surface, 

where the filter extends outside the fault surface. In our case, the fault surface is cut out from a 

much larger stress distribution. So stress values outside the fault surface are available, and we 

use them when applying the filter near the edges of the fault surface.

Scaling

To perform dynamic rupture simulations, the stresses must be scaled to be consistent with 

the friction law. In the high-stress (center) portion of the fault surface, the stresses must be high 

enough to sustain the rupture but belo the fault’s yield stress. In the lo -stress (side and 

bottom) portion of the fault surface, the stresses must be low enough to stop the rupture.

Our recommended procedure is to generate a group of 10 or 20 stress distributions, and 

then scale all of them at once. Then the scaling formula is:

 

    (289)

where is the nth stress field, and the numbers and are selected as needed to make the 

stress field consistent with the friction law. Notice that the same scale factor is used for all the 

stress distributions in the group; this allows the amplitude of the stress fluctuations to vary 

among the stress realizations.

After scaling, the stress field may exceed the fault’s yield stress at some points. If this 

occurs, then at those points either the stress field should be reduced, or the friction parameters 

should be modified, to increase the yield stress.

See the example at the end of this appendix for further discussion of scaling.
87

Hypocenter Selection

When heterogeneous initial stresses are used, the hypocenter should be located within a 

region of relatively high initial stress.

In our computer code, the user specifies a range of permissible hypocenter locations, 

along-strike and along-dip. For each permissible location , the code calculates a weighted 

average of the initial stresses at nearby points, like this:

∑ ∑ 

  

 

 

 

 

(290)

where and are the spacing between adjacent points on the fault surface, and where is the 

shear stress (after scaling). The coefficients are called the hypocenter selection mask. 

We typically take them to be Gaussian, which is defined as:

 

 ( 

 

 

  

 

 

 

 ) (291)

The constant is chosen so that

∑ ∑ 

 

 

 

 

 

 (292)

Typically we take and . However, the precise form of the 

hypocenter selection mask is not important.

The location on the fault with the highest weighted-average nearby stress is 

chosen as the hypocenter.

Size of Nucleation Zone

The main report describes a new nucleation method, which combines a small zone of 

forced rupture (equation 77) and a larger zone of reduced fracture energy (equation 78). Both 

zones are circles centered at the hypocenter. If this method is to be used, then it is necessary to 

select the radius of the zone of reduced fracture energy. (The radius of the zone of forced 

rupture is then .)

We choose the radius of the zone of reduced fracture energy to be

 F NU (293)

where is the Day radius (equation 76), and is a user-specified parameter called the 

nucleation factor. The nucleation factor is chosen to be large enough, so that the nucleation 

method operates reliably. In our tests, we find that values in the range 1.00 to 1.25 work 

well.

Some experimentation may be necessary to select an appropriate value of , by 

running dynamic rupture simulations with different values of . If is too small, the 
88

simulated rupture will die out before exiting the nucleation zone, or shortly after exiting the 

nucleation zone. If is just slightly too small, the simulated rupture may slow down or stop 

near the edge of the nucleation zone, then speed up and continue after a delay. When is 

large enough, then the rupture should propagate smoothly through the nucleation zone and on 

into the surrounding fault. On the other hand, should not be excessively large, because we 

desire to minimize the effect of the nucleation process on the rupture propagation. By performing 

several dynamic rupture simulations, with several different values of , and observing the 

resulting simulated rupture propagation, it is possible to select an appropriate value of .

The Day radius (Day, 1982) is defined to be

 

 

 

 

 

 

 (294)

Here, is the shear modulus, is the initial shear stress, is the static yield stress, is the 

dynamic sliding stress, and is the slip-weakening critical distance, that is, the amount of slip 

over which the frictional stress declines from to . The fracture energy per unit area is 

  .

To compute the Day radius, we use the approximation

 

 

 

 ( avg )

 

 

 (295)

Here, is the shear modulus in the neighborhood of the hypocenter, is the static coefficient of 

friction at the hypocenter, is the dynamic coefficient of friction at the hypocenter, is the 

slip-weakening critical distance at the hypocenter (not taking into account the reduced critical 

distance used for nucleation), and is the initial normal stress at the hypocenter. The fraction 

 is the average ratio of initial shear stress to initial normal stress, over a user-specified 

rectangular portion of the fault surface. Because of the scaling by normal stress, which is 

described below, is actually proportional to the “initial stress” value at this stage of the 

post-processing.

Note that this calculation implicitly assumes the use of a slip-weakening friction law. 

Such a law is characterized by a static coefficient of friction , a dynamic coefficient of friction 

 , and a slip-weakening critical distance , which is defined to be the amount of fault slip over 

which the coefficient of friction decreases from to .

Depth Conditioning

If the goal is to model earthquakes large enough to rupture the full depth of the 

seismogenic zone, then depth conditioning is used to limit rupture extent at the bottom of the 

fault. This is appropriate because the base of the seismogenic zone is defined by a systematic 

change in the properties of the fault, such as a transition from brittle to ductile behavior.

Depth conditioning systematically reduces the shear stress with increasing depth. The 

shear stress is multiplied by a depth-conditioning function, , which varies as a function of depth

(equation 69):

 

   (296)
89

In the examples shown here, the depth conditioning function is taken to be piecewise 

linear (equation 70):

 

 {

 for m

 m for m m

 m for m m

(297)

It is of course possible to make the depth conditioning function be smooth rather than piecewise 

linear.

If the goal is to model smaller earthquakes that do not rupture the full depth of the 

seismogenic zone, then depth conditioning should not be used. In this case, the selection 

template is used to control rupture extent both vertically and horizontally.

Scaling by Normal Stress

A dynamic rupture simulation will assume some pattern of initial normal stress on the 

fault. In general, friction forces are proportional to the normal stress. So, to remain consistent 

with the friction law, we scale the shear stress by a factor proportional to the normal stress:

 

   (298)

where is the initial normal stress on the fault at point , and where is some 

constant chosen to attain consistency with the friction law.

Another way of saying this is that we treat our randomly generated values as if they lie in 

“friction coefficient space.”

Example of Scaling

We present an extended example of scaling, to illustrate the group scaling (equation 289) 

and normal scaling (equation 298) procedures. The scaling technique described here was used in 

the examples that accompany this report. We emphasize that this is an example, that there are 

other acceptable ways to perform the scaling, and that, in any event, the scaling must be adjusted 

to produce stress fields consistent with the chosen friction parameters.

1. For this example, we assume a slip-weakening friction law, with static coefficient of 

friction and dynamic coefficient of friction 

2. For this example, we assume the initial normal stress has a constant value 

 

 .

3. Because of the normal stress scaling, the value of 

  appearing on the left side 

of equation 298 should be regarded as lying in “friction coefficient space,” that is, it 

is the ratio of the initial shear stress divided by the normal stress (and hence its values 

should lie in the same range as the friction coefficients).

4. Recall that the selection template (equation 73) specifies high-stress and 

low-stress regions of the fault surface. In figure 8, the high-stress region is the deepest 

red area, and the low-stress region is the deepest blue area.
90

5. Let ̅ 

 denote the average value of 

  in the high-stress region. Let ̅

 

denote the average value of 

  in the low-stress region.

6. For this example, we select a target value for ̅ 

 

, which is 

 . 

However, we do wish to allow some random variation in the value of ̅ 

 , that is, we 

do not want to force ̅ 

 to be exactly equal to  

.

If we define the earthquake strength ratio (Day, 1982) to be ( ̅ 

 ) 

( ̅ 

 ), then the target value of is (  ) (  ) (this is 

how we selected the numerical value of 

 ). The earthquake strength ratio is an 

indicator of how fast the rupture propagates, with lower ratios yielding faster 

propagation. By allowing some random variation in value of ̅ 

 , we are thereby 

allowing random variation in and hence in the rupture propagation speed. The 

reason that we employ group scaling (equation 289) is to obtain this random variation 

in rupture speed.

7. For this example, we select a target value for ̅

 , which is 

 . The value 

must be chosen far enough below so that the rupture comes to a stop when it 

enters the low-stress region. For this example, we have chosen the numerical value of 

 

 to satisfy 

 

 . Again, we want to allow some random 

variation in the value of ̅

 .

8. The parameter in equation 289 is chosen so that the average value of ̅ 

 ̅

 is 

made equal to 

 

 .

9. Then, each parameter in equation 289 is chosen so that 

 

( ̅ 

 ̅

 ) .

10. Then, equation 289 is used to calculate all the values of 

  .

11. Then, we adjust the values of 

  so they are not excessively high or low. In 

this example, we impose the constraint 

  . Any computed 

value of 

  that lies outside this range is adjusted to bring it into range.

12. Then, for points within the nucleation zone, adjust the values of 

  so they are 

above the dynamic coefficient of friction. This is done so that there are no low 

stresses within the nucleation zone that could interfere with the nucleation process. In 

this example, we impose the constraint 

  for points within the 

nucleation zone. Any computed value of 

  that lies outside this range is 

adjusted to bring it into range. The value was obtained as being approximately 

 ( 

 ).

13. Then, we increase the static coefficient of friction, where necessary, to make the yield 

stress higher than the initial stress. In this example, we impose the constraint 

 

  . At any point where this constraint is violated, we increase 

so 

the constraint is satisfied (thereby producing a non-uniform value of on the fault 

surface).

14. Finally, we compute the actual initial shear stresses by applying equation 298 with the 

parameter equal to .
91

The end result of this scaling procedure is a heterogeneous initial stress field that is 

consistent with the chosen friction law. In this example, the initial stress in the high-stress region 

is consistent with an average strength ratio ; the initial stress in the low-stress region is 

low enough to arrest the rupture; there is random variability in the strength ratio; the initial stress 

is everywhere below the yield stress; and the initial stress inside the nucleation zone is above the 

dynamic sliding stress.
92

Appendix H. One-Point Statistics

In this appendix we present some empirical results on the one-point statistics produced by 

our method and the effect of applying a low-pass spatial filter to the stress field.

The term one-point statistics refers to the statistical properties of the stress field evaluated 

at individual points on the fault surface (such as, not taking into account any correlations 

between stress values at different points on the fault surface). It has been suggested (for example, 

Lavallée and others, 2006) on the basis of slip inversions that the one-point statistics should be 

non-Gaussian and, in particular, that it should be more heavy-tailed than a Gaussian. We present 

evidence that our method produces heavy-tailed one-point statistics.

To assess the one-point statistics, we treat the all the stresses as being values of a random 

variable and perform a simple test to see if their distribution appears Gaussian or heavy-tailed. 

To perform the test, we partition the range of stress values into 50 bins. The fraction of stress 

values that fall into each bin is taken to be an estimate of the probability density function. Then, 

we find both the best-fitting Gaussian distribution and the best-fitting Cauchy distribution. Note 

that the Cauchy distribution is so heavy-tailed that its variance is infinite. Then we calculate the 

ratio , where is the residual from the Cauchy fit, and is the residual from the Gaussian 

fit. With this metric, high values tend to indicate a Gaussian distribution, and low values tend to 

indicate a heavy-tailed distribution. We performed this calculation for 20 random stress 

distributions.

Recall that in our method, the stress pattern on the fault surface is cut out from a large 

stress distribution that is thousands of times larger than the fault surface. If we use stress values 

from the entire large stress distribution, before the fault surface is cut out, the average ratio is 

38.5. This indicates that the one-point statistics of the large stress distribution is very well 

described as being Gaussian.

If we use stress values from just the selected fault surface, the average ratio is 1.26 if the 

stress is not low-pass filtered, or 1.23 if the stress is filtered. If we use just the stress values from 

a rectangular portion of the fault surface that lies in the high-stress area, the average ratio is 1.90 

if the stress is not filtered, or 1.85 if the stress is filtered. In either case, the ratios indicate that the 

one-point statistics are more heavy-tailed than a Gaussian. It is interesting to note that applying a 

low-pass spatial filter does not alter the heavy-tailed nature of the one-point statistics.

Discussion

Our method, like some prior methods (for example, Andrews and Barall, 2011), generates 

stress values by summing large numbers of random values. Such methods have a strong tendency 

to generate Gaussian one-point statistics, due to the action of the Central Limit Theorem. This is 

illustrated by the nearly perfect Gaussian one-point statistics that are observed on the large stress 

distribution.

There has been some question as to whether it is possible for methods of this type to 

produce heavy-tailed, non-Gaussian statistics. Here we present evidence that this is indeed 

possible, by selecting just a portion of the generated stress values. Although the entire set of 

stress values on the large stress distribution is always Gaussian, the stress values on the selected 

fault surface can have non-Gaussian statistics.
93

Table 1. Mathematical notations.

The following table is a list of mathematical notations used in the main portion of this

report. (Note that some appendixes may use slightly different notations.)

Notation Definition

 Total area of all stress perturbations with radius . See equation 47.

 

 Filter kernel, used during post-processing to apply a low-pass filter to the 

stress field. See equation 286 in appendix G.

 Slip-weakening critical distance, that is, the amount of slip over which the 

friction drops from its static value to its dynamic value.

 

 Kernel, used during post-processing to select a hypocenter location with 

high local stress. See equation 290 in appendix G.

 

 Elastic energy of the basis displacement field for wave vector . 

See equation 15.

 

 Maximum allowed elastic energy of the basis displacement field for wave 

vector ; when applying state restriction. See equation 54.

 Total elastic energy of a linear combination of the basis displacement 

fields. See equation 16.

 

 Elastic energy density of the basis displacement field for wave vector 

 . Implicitly a function of position in three-dimensional space. See 

equation 130 in appendix B.

 Total elastic energy density for a linear combination of the basis 

displacement fields. Implicitly a function of position in three-dimensional

space. See equation 147 in appendix B.

 

 Amplitude of basis displacement field for wave vector . See 

equations 8–10.

 

 Maximum allowed amplitude of the basis displacement field for wave 

vector ; when applying state restriction. See equation 53.

 Nucleation factor, which is a multiplier applied to the Day radius to obtain 

the radius of the zone of reduced fracture energy. Used in the two-stage 

nucleation method. See equation 293 in appendix G.

 Energy rolloff function. Gives the expected value of the energy 

 

of a basis displacement field, as a function of the scalar wavenumber . 

See equations 31, 49, and 51.

 Variable used in calculating 

 . See equation 56.

  Wave vector in the fault plane, associated with the pattern of fault slip.

 Wave number in the vertical (along-dip) direction along the fault plane.
94

Notation Definition

 Wave number in the horizontal (along-strike) direction along the fault 

plane.

 Scalar wave number, 

 

  

.

 

 Nyquist wave number in the vertical (along-dip) direction. See equation 

62.

 

 Nyquist wave number in the horizontal (along-strike) direction. See 

equation 63.

 Fault size vertically. Also the vertical period of the periodic boundary 

conditions.

 Fault size horizontally. Also the horizontal period of the periodic 

boundary conditions.

 

Integer identifying a grid point in the vertical direction. See equation 58.

 

Integer identifying a grid point in the horizontal (along-strike) direction. 

See equation 59.

 

Integer representing a number of periods in the vertical direction. See 

equations 17 and 60.

 

Integer representing a number of periods in the horizontal (along-strike) 

direction. See equations 17 and 61.

 

Integer representing a number of grid points in the vertical direction.

 

Integer representing a number of grid points in the horizontal (alongstrike) direction.

 Rate of occurrence of earthquakes with radius between and . See 

equation 42.

 Number of stress perturbations with radius between and that 

exist at an instant of time. See equation 46.

 Selection template. Used in the selection algorithm, it identifies the 

portions of the fault surface where high or low stress is desired. 

 Day radius, the minimum radius a rupture must attain so that further 

expansion of the rupture is energetically favorable. See equation 76. Also 

see equations 294 and 295 in appendix G.

 Radius of the zone of forced rupture, used in the two-stage nucleation 

method. See equation 77.

 Radius of the zone of reduced fracture energy, used in the two-stage 

nucleation method. See equation 78. Also see equation 293 in appendix G.

 

 Amplitude of generalized basis displacement field for wave vector 

 . See equations 82–84 in appendix B.
95

Notation Definition

 

 Rotation template, used during post-processing to decide whether or not to 

rotate the stress field by 180 degrees. See equation 285 in appendix G.

 

  Basis displacement field, which is an approximate solution to the 

elastostatic fault slip problem with sinusoidal fault slip. Displacement in 

coordinate direction as a function of position in space. There is a 

separate basis displacement field for each admissible wave vector, hence 

 is also implicitly a function of wave vector . See equations 8–

10.

 

 Total displacement field, obtained by forming a linear combination of the 

basis displacement fields. Implicitly a function of position in threedimensional space. See equation 146 in appendix B.

 Random variable, uniformly distributed between 0 and 1.

 Random variable, uniformly distributed between 0 and 1.

 Coordinates in three-dimensional space.

 Vertical coordinate.

 Horizontal coordinate along-strike.

 Horizontal coordinate perpendicular to the fault plane.

 Real-valued random variable that has the expected one-point statistics of 

the filtered stress field . See equation 187 in appendix 

D.

 Plane wave amplitude for the shear stress generated by the basis 

displacement field for wave vector . See equation 182 in appendix 

D.

 Decay exponent of generalized basis displacement field. Implicitly a 

function of position perpendicular to the fault plane and scalar 

wavenumber . See equation 90 in appendix B.

 In thermodynamics, a parameter inversely proportional to temperature. In 

fault mechanics, a parameter that determines the correlation length and the 

overall level of stress fluctuations.

 Strain tensor. Implicitly a function of position in three-dimensional space.

 Parameter appearing in the energy rolloff function , adjusted to set 

the correlation length. See equations 49 and 51.

 Parameter appearing in the energy rolloff function . See equation 51.

 Complex phase of the convolution 

 

 ̃ ̃  at the point where 

its absolute value is maximized. See equation 74.
96

Notation Definition

 Correlation length. Defined to be ̂, where ̂ is chosen to 

maximize the value of ( ̂) ̂ .

 Volumetric Lame parameter.

 Shear modulus. Also the shear Lame parameter.

 Phase of basis displacement field for wave vector . See equations 

8–10.

 Total phase of generalized basis displacement field. Implicitly a function 

of position on the fault plane and wave vector . See 

equation 89 in appendix B.

 Stress tensor. Implicitly a function of position in three-dimensional space.

 Normal stress acting on the fault plane. Implicitly a function of position on 

the two-dimensional fault plane.

 Shear stress acting on the fault plane. Implicitly a function of position on 

the two-dimensional fault plane. We take the shear stress to be oriented 

along-strike.

 

 Shear stress acting on the fault plane, due to the basis displacement field 

for wave vector . It is implicitly a function of position on 

the fault plane. See equation 20.

 

 

  Mean-square value of the shear stress acting on the fault plane, due to the 

basis displacement field for wave vector . See equation 21.

 

 Total shear stress acting on the fault plane, due to a linear combination of 

the basis displacement fields. See equation 68.

 

 Mean-square value of the shear stress acting on the fault plane, due to a 

linear combination of the basis displacement fields. See equation 22.

 

 

  Complex-valued total shear stress acting on the fault plane, due to a linear 

combination of the basis displacement fields. Used in the selection 

algorithm. See equation 72.

 Total shear stress acting on the fault plane, due to a linear combination of 

the basis displacement fields, for wavenumbers . See equation 185 in 

appendix D.

 

 

 Mean-square value of the shear stress acting on the fault plane, due to a 

linear combination of the basis displacement fields, for wavenumbers 

 . See equation 186 in appendix D.

 Depth conditioning function, which systematically reduces the shear stress 

with increasing depth. See equation 69. Also see equation 296 in appendix 

G.
97

Table 2. Summary of thermodynamic analogies.

The following table summarizes our analogy between thermodynamics and earthquake 

mechanics. Each row in the table lists a concept from thermodynamics and the analogous 

concept for earthquake faults.

Thermodynamic concept Earthquake fault concept

Physical system. A planar fault, with periodic boundary 

conditions along-strike and along-dip, and the 

surrounding elastic medium.

Micro-state—A complete specification of the 

system’s microscopic state.

A pattern of slip on the fault surface and the 

resulting distortion of the surrounding elastic 

medium.

Energy of a micro-state. Total elastic energy due to the distortion of 

the surrounding elastic medium.

Degrees of freedom—Partitioning a system 

into a set of subsystems, so that the energy of 

the complete system equals the sum of the 

energies of the subsystems. Each subsystem is 

a degree of freedom.

Orthogonal basis displacement fields—

Decomposing an arbitrary pattern of slip on 

the fault into a sum of basis displacement 

fields, so that the total elastic energy of the 

complete slip pattern is equal to the sum of 

the elastic energies of the basis displacement 

fields.

Density of states—Number of micro-states 

per unit volume of phase space, or per unit 

energy interval.

Density of states—Number of basis 

displacement fields with wave vectors that 

satisfy periodic boundary conditions, per unit 

volume of wave-vector space.

Equipartition theorem—Each degree of 

freedom has the same expected energy (under 

certain conditions).

Equipartition theorem—Each basis 

displacement field has the same expected 

elastic energy (under certain conditions).

Ultraviolet divergence—A mathematical 

divergence that arises in systems with 

infinitely many short-wavelength degrees of 

freedom, because equipartition then predicts 

infinite short-wavelength energy. For 

example, blackbody radiation is predicted to 

contain an infinite amount of ultraviolet 

radiation.

Ultraviolet divergence—A mathematical 

divergence that arises because there are 

infinitely many short-wavelength basis 

displacement fields, and equipartition then 

predicts an infinite amount of shortwavelength elastic energy and stress.
98

Thermodynamic concept Earthquake fault concept

Quantization—Reduces the number of microstates by requiring that some physical 

quantities (such as energy) can only assume 

discrete values instead of being continuously 

variable. Can resolve an ultraviolet 

divergence by reducing the energy in shortwavelength degrees of freedom to below the 

equipartition value.

Finite fault strength—Reduces the number of 

acceptable patterns of slip on the fault surface 

by requiring that the shear stress cannot 

exceed the yield strength of the fault. (We 

enforce this rule in the weak, or mean-square, 

sense; rather than enforcing it at each point on 

the fault surface.) Resolves the ultraviolet 

divergence by reducing the energy in shortwavelength basis displacement fields.

Planck distribution—For blackbody radiation, 

a formula that specifies how the expected 

energy for each degree of freedom falls off 

toward zero as the wavenumber increases.

Energy rolloff function—A formula that 

specifies how the expected energy for each 

basis displacement field falls off toward zero 

as the wavenumber increases.
99

Figures. 

Figure 1A.

Dynamic 

rupture 

code 

Fault 

slip

Seismic 

waves 

Initial 

stress 

conditions 

Fault 

geometry 

Material 

properties 

Friction 

law
100

Figure 1B.

Kinematic

rupture 

code 

Seismic 

waves 

Fault 

slip

Fault 

geometry 

Material 

properties 

Figure 1. Components of dynamic and kinematic earthquake simulations. (A) A dynamic rupture 

simulation models the frictional processes that occur on the fault, so the inputs include the friction law 

and initial stresses, and the simulation code must calculate both fault slip and seismic waves. Our 

method is a way to inject randomness into the initial stress conditions. (B) A kinematic rupture simulation 

treats the fault slip as being known a priori, so there is no need to specify a friction law or initial stresses, 

and the simulation code only needs to calculate the seismic waves.
101

Figure 2. Two-dimensional Ising model. The Ising model is a model of a ferromagnetic material such as 

iron. The circles represent atoms in a crystal lattice, each of which can interact only with its four nearest 

neighbors. Each atom has a “spin,” shown by an arrow, that can point either “up” or “down.” The spin 

can be thought of as being a small magnet. The interaction energy between two nearest-neighbor atoms 

depends on whether their spins are parallel or antiparallel. The Ising model illustrates how a system can 

self-organize into large structures such as magnetic domains. 
102

x 

x 

x 

L 

L 

Figure 3. Three-dimensional coordinate system. We assume a vertical strike-slip fault. The fault lies in 

the -plane, with the axis pointing up and the axis pointing along-strike. The axis is 

perpendicular to the fault. We impose periodic boundary conditions in the and directions (but not 

the direction), with periods and , respectively. So the area of the fault is . 
103

Viewport 

T 

Figure 4. Blackbody radiation. A container is maintained at a temperature , and a small viewport 

allows the experimenter to observe the electromagnetic radiation inside the container. The problem is to 

compute the spectrum and energy density as a function of temperature. The equipartition theorem 

predicts infinite energy density, a difficulty known as an ultraviolet divergence. However, quantum 

mechanics greatly reduces the number of possible states of the system, with the result that the energy 

density is rendered finite. 
104

Figure 5A.
105

Figure 5B 

Figure 5. Energy rolloff function. (A) Energy rolloff function plotted as a function of the scalar 

wavenumber . The plot begins with and decreases monotonically, approaching in the limit 

 . As explained in the main text, the finite strength of the fault restricts the possible states of the 

system, and, as a result, the expected energy of a displacement field with wavenumber is given by 

 where is the thermodynamic parameter analogous to reciprocal temperature. By 

comparison, equipartition predicts that the expected energy is independent of , which leads to an 

ultraviolet divergence. (B) One-dimensional power spectral density function 

 

 plotted as a 

function of . The area under this curve must be finite, for our method to yield a finite level of stress. The 

maximum value occurs at where is defined to be the correlation length. The plots in this 

figure use 

 

 

 where . 
106

Figure 6A.
107

Figure 6B.

Figure 6. Examples of random stress distributions generated by our thermodynamic method. Levels of 

stress on the fault surface are shown by varying colors, which indicate the ratio of shear stress divided 

by normal stress. Rupture extent is limited by the low-stress areas along the sides and bottom of the 

fault. (A) Stress pattern used in Southern California Earthquake Center Dynamic Rupture Code 

Validation benchmark TPV16. (B) Stress pattern used in benchmark TPV17. 
108

Stress distribution: 

2,457.6 km square, 

32,768 × 32,768 

nodes. 

Fault surface: 

48.0 km × 19.5 km, 

641 × 261 nodes. 

Figure 7. Selection process for controlling rupture extent along-strike. A random stress distribution is 

generated over an area thousands of times larger than the fault surface. Then, the entire large area is 

searched to find a portion, the size of the fault surface, which by chance has high stress in the center 

and low stress at the sides. In this example, the large area is 32,768 grid points square, which is 2,457.6 

kilometers (km) square assuming a grid point spacing of 75 meters (m). The fault area is 641 by 261 

grid points, or 48.0 km by 19.5 km. The figure is not to scale. If drawn to scale, the fault surface would 

have about one-tenth the area shown.
109

Figure 8. Selection template. Colors show the value of the selection template , which is used 

in the selection process. It is positive in the area where high stress is desired, in this example the central 

30 kilometers (km) of the fault surface. It is negative where low stress is desired, which is to the left and 

right sides. It is zero where the stress level is unimportant, in this example the bottom 2.5 km of the fault 

surface where rupture extent is controlled by depth conditioning. There is a gentle taper between the 

positive and negative areas, spanning 6 km, to allow for an irregular border between high and low stress 

areas. The selection template tapers to zero over a span of 0.5 km, part of which lies outside the fault 

surface (not shown). 
110

Figure 9A.

Figure 9B.
111

Figure 9C.

Figure 9. Effect of varying the correlation length. Colors show the level of stress on the fault surface, 

computed for three different correlation lengths. Images are shown without the use of depth 

conditioning. (A) Correlation length = 10 kilometers (km). This is too small, as shown by the fact that it is 

impossible to find a high-stress area 30 km across. (B) Correlation length = 100 km. This is too big. 

Notice that the borders between the high- and low-stress areas are fairly well-defined. (C) Correlation 

length = 30 km. This is the right size, in this example where we seek a high-stress area 30 km across. 

Note the ragged borders between the high- and low-stress areas. Our solution is to choose the smallest 

correlation length that gives a sufficiently large area of high stress to sustain the simulated earthquake 

rupture. 
112

Figure A–1. Slider-block model of an earthquake fault. The fault is modeled as a rectangular array of 

blocks, each connected by springs to its four nearest neighbors. The blocks can slide within the plane, 

subject to the forces produced by the springs, and the displacement of a block from its initial position 

represents fault slip. Each block is also subject to a tectonic driving force, which could be supplied by a 

fifth spring, which is connected to a rigid plate (not shown in the figure). Although the slider-block model 

has an obvious similarity to the Ising model, it is not suitable for our purposes because it has incorrect 

energy. If a sinusoidal displacement pattern is applied to the blocks, the energy stored in the springs is 

proportional to 

 

, where is the wavenumber, whereas sinusoidal slip on a fault produces elastic 

energy proportional to .