1

A RECIPE TO GENERATE PHYSICALLY CONSISTENT EARTHQUAKERUPTURE MODELS FOR GROUND MOTION PREDICTION
Mariagiovanna Guatteri ,1 ,P. Martin Mai 2 ,Gregory C. Beroza 3
1

Swiss Reinsurance America, Armonk, NY 10504
Institute of Geophysics, ETH Hoenggerberg, Zurich, CH-8093, Switzerland
3
Department of Geophysics, Stanford, CA 94305-2215
2

Summary
Accurate predictions of the intensity and variability of ground motions from future large
earthquakes depend strongly on our ability to simulate realistic models of the earthquake source.
In order to make these simulations more accurate, we have developed a procedure to generate
physically consistent earthquake rupture models. We term these models pseudo-dynamic because
they are kinematic models that are designed to emulate important characteristics of fullly
dynamic rupture models. The user-specified parameters to generate pseudo-dynamic rupture
models are the earthquake magnitude and hypocenter, with the option of specifying the fault
dimensions as well. We construct pseudo-dynamic models by first generating a slip distribution
as a realization of a spatial random field that is consistent in its overall scaling and spatial
variability with slip distributions observed in past earthquakes[Mai and Beroza, 2000; Mai andBeroza, 2002]. We compute the static stress drop associated with this slip distribution, which inturn is used to assign the temporal evolution of slip (rupture velocity, rise time, slip-velocity
function) by applying a series of empirical relationships derived from the analysis of a set of
spontaneous rupture models. A simple energy-budget calculation is used to discard source models
that are not realizable as spontaneous ruptures.
Our approach is based on the findings of Guatteri et al. [2003], incorporating ideas developedby Mai et al. [2001] and Guatteri et al. [2002] for a source characterization that adequately
describe the dynamics of rupture without having to do full dynamic simulations. The pseudodynamic approach circumvents the limits imposed by the extreme computational demand ofgenerating fully dynamic rupture models for simulating multiple realizations of a scenario
earthquake for strong-ground motion prediction. While the relationships between source
parameters described in this paper are significant simplifications of the true complexity of thephysics of the rupture process, they help identify important interactions between source properties
that are relevant for strong ground motion prediction.

Introduction
Statistical models have long been proposed to characterize heterogeneity in earthquake slip andto relate it to observable characteristics of earthquake ground motion [Haskell, 1966; Aki, 1967].
Andrews [1980b] presented a theoretical formulation for this problem for heterogeneous faults
that has formed the foundation for most subsequent efforts [Herrero and Bernard, 1994; Frankel,
1991; Zeng et al., 1994]. Andrews [1980b] showed that a two-dimensional slip function D(x,z)
with a fractal dimension D = 2 results in a far-field spectral decay of displacement as ω-2. The
fractal dimension D is related to the wavenumber decay of the two-dimensional Fourier spectrumof the slip function, D(k) (where k is the wavenumber); for D = 2, the wavenumber spectra
decays as k-2. Using this relation, Herrero & Bernard [1994] propose the “k-square” model that
introduces a source-size dependent length scale Lc as well as a wavenumber-dependent behaviorof the rise time distribution. Zeng et al. [1994] develop the composite source model based on themodel of Frankel [1991] in which elementary circular sources with a fractal size distribution aresummed to form the complete two-dimensional slip function. An explicit assumption in these2models is that stress drop is scale-independent. Moreover, none of these methods directly takes
into account slip heterogeneity as imaged using seismic data in finite-source models.
In order to characterize the slip complexity of past earthquakes, Somerville et al. [1999] adopt a
deterministic approach to correlate the size and number of asperities with seismic moment for aset of finite-source rupture models. Mai and Beroza [2002], in contrast, use a spatial random-fieldmodel to describe the complexity of earthquake slip distributions. They find that a von Karman
autocorrelation function for which the correlation length increases with source dimension mostclosely represents the spectral properties of existing earthquake slip models. Increasing
magnitude and hence increasing source dimensions with constant correlation lengths would meanthat the rupture would be comprised of many isolated high-slip (high stress-drop) asperities,
which would show extremely large slip (and stress drop) in order to accommodate the seismic
moment. Unless the correlation length grows with earthquake size, these areas of positive stressdrop will be too weakly connected for the rupture to propagate spontaneously because they will
be separated by large areas of negative stress drop, and therefore negative elastostatic energy,
which will act to stop the rupture.
Guatteri et al. [2002] find that including spatial and temporal variations in slip, slip rise time,and rupture propagation that are consistent with dynamic rupture models exerts a strong influenceon near-source ground motion and has the potential to improve strong ground motion prediction.
Their results lead to a feasible approach to specify the variability in the rupture time distribution
in kinematic models consistently with dynamic source properties. Mai et al. [2001] and Guatteri
et al. [2002] calculate ground motion from a “first-generation” set of pseudo-dynamic models,demonstrating the potential of such a procedure. In this paper, we focus on the generation of thesource models by expanding and improving upon the previous “first-generation” pseudo-dynamic
source characterization, and leave the ground motion calculation for future applications.

Pseudo-Dynamic Source Characterization
In this paper we define as a pseudo-dynamic source model a kinematic model in which therelevant source parameters (slip, rupture velocity and slip-velocity function) are specified in sucha way that they emulate both slip distributions of past earthquakes and the temporal behavior of
spontaneously propagating dynamic rupture models. We developed the pseudo-dynamic approach(henceforth abbreviated as PD) by investigating the relationships between kinematic and dynamic
source parameters for a set of dynamic rupture models representing strike-slip earthquakes with amagnitude range of 6.4 < M < 7.2.
The complexity of the spatial slip distribution in our PD models is characterized as a spatial
random field [Mai and Beroza, 2002] that is stochastically consistent with slip models for pastearthquakes. The starting point of the PD model consists of generating a spatial distribution ofslip as a realization of a spatial random field.
The characterization of the temporal slip evolution in finite-source rupture models for pastearthquakes is limited by a lack of resolution of the variability of temporal slip parameters, suchas rupture velocity and slip rise time. The fundamental idea of our approach is that we canovercome this limitation by assuming that dynamic rupture modeling can provide physicalconstraints to the temporal slip evolution characterization, given a spatial slip distribution. We
therefore perform spontaneous dynamic rupture simulations to develop the relationships between
kinematic and dynamic rupture parameters to constrain the temporal rupture characteristics. The
following section describes our dynamic modeling approach and the various steps towarddefining physically consistent rupture velocity, rise time and slip-velocity functions.

Initial Dynamic Rupture Modeling
The first PD source characterization was developed by Mai et al. [2001] and improved by
Guatteri et al. [2002]. Both their results were based on the set of 6 stochastic-dynamic models

3developed by Guatteri et al. [2003] computed using a Boundary Integral Method [Boatwright andQuin, 1986; Das and Kostrov, 1987; Quin and Das, 1989]. These models were discretized into agrid size of 0.75 km, allowing for a maximum frequency of about 2 Hz (for this method), and aslip-weakening fault constitutive relationship was assumed [Andrews, 1976; Day, 1982] (Figure1). Table 1 summarizes and defines the relevant dynamic parameters together with a brief
explanation of how we assign or calculate them to generate our set of stochastic-dynamic models.

Figure 1: Slip-weakening model showing the relevant energy contributions: elastostatic energy Eel (left),fracture energy Efr and relaxation work Erx (center). The right panel displays the case of negative stressdrop Δτ for which fracture energy is still defined. See Table 1 for further explanations. τ1 is the final
stress, which may be different from the dynamic frictional level.

TABLE 1: Relevant dynamic parameters
Parameter
1.
σnGeneral Information
Description
Normal stressValue
150 MPa

In this study
Comment
Uniform over the fault

2.

µdCoefficient of dynamic friction

0.43

Uniform over the fault

3.

τf

µdσn

Uniform over the fault

4.

Δτ

Dynamic resistance(frictional stress level)
Stress drop

τ0 - τtCalculated from slip
distribution (eq. 3)

5.
6.

τ0
Dc

Initial shear stressSlip-weakening distanceΔτ + τf
Sampled from [0.25 – 0.4] m forthe whole fault

7.

Efr

Apparent fracture energy(also termed Gc)½ Dc (τy - τf )8.

τy- τ0Strength excess τex

Calculated from Efr and Dc values

9.
10.

S
τy

(τy - τ0 ) / Δτ
µsσn

11.
12.
13.

µs
Eel
Ers

Strength parameter
Static resistance(upper yield stress)
Coefficient of static friction
Elastostatic energy
Radiated seismic energyτs / σn- ½ sliptot (τf - τ0 )Eel – Efr;
Erx = 0 if the frictional stress τf isequal to the final stress τ1

Assigned to lower boundwhere Δτ < 0Simulated from stressdrop distribution using(eq. 6)Lower bound τex = 0.5
where Δτ <0. Upper
bound given by Smax
Upper bound Smax = 2.5
Calculated from (eq. 8)
and (eq. 5)

Assumption: static stressdrop is equal to the
dynamic stress drop

4
We assume the same velocity profile as in Guatteri et al. [2003] derived from Boore andJoyner [1997]. We force the rupture to nucleate at a given hypocenter location and to propagate ata constant rupture speed within a nucleation area whose size is calculated according to Day[1982].
Our procedure to generate the new set of stochastic-dynamic models is very similar to that of
Guatteri et al. [2003]. We first generate our target slip distribution using the method of Mai andBeroza [2002], from which we compute the corresponding static stress drop using the method ofAndrews [1980]. Once we set the hypocenter location, instead of assigning the fracture energy
distribution by the trial and error procedure as in Guatteri et al. [2003], we applied the empirical
relationship between fracture energy and stress drop and crack length developed by Guatteri et al.
[2002]. In a later section in this paper we describe this relationship in detail and discuss how we
modified it to satisfy the assumption of sub-shear average rupture velocity over the entire fault
plane for the stochastic-dynamic models.
The set of stochastic-dynamic rupture models are representative of strike-slip earthquakes
with a magnitude range of 6.4 < M < 7.2. Table 1 summarizes the important features of the
dynamic parameters, while Figure 2 schematically shows the roadmap of our PD proceduredevelopment. Table 2 summarizes the parameters used to define a PD model.
Relationship between
rupture velocity, v,
Gc, Δτ, and Lh
Equations (4,5)

Second-Generation pseudo-dynamic modelingFirst-Generation pseudo-dynamic modeling

Stochastic-Dynamic models[Guatteri et al., 2003]

Δτ from slip

Mai et al., [2001]•
Gc = constant•
Rise time, τr, based onhealing front propagating at
shear-wave velocity

Guatteri et al., [2002]•
Gc simulated, Equation (6)•τr simulated, Equation (10)Slip-velocity function parameterization as aquasi-dynamic Kostrov function

THIS STUDY

Stochastic-dynamic simulation using Gc from Equation (6)••New parameterization of τr
New parameterization of SVF, including peak slip-velocity, Vmax and time pulse, Tp

5

Figure 2: (previous page) Roadmap showing the various phases involved in the development of the pseudo-dynamic
procedure presented in this paper.

TABLE 2: Kinematic, dynamic and pseudo-dynamic parameters
Slip
Lh, distance from the hypocenter

Rupture characteristics

Assigned kinematic parametersΔτ, stress dropGc, fracture energy

Deduced from kinematic, but not
rupture characteristics

Intermediate, dynamic parametersv, rupture velocity
τr, rise time
SVF (Vmax, Tp)

Rupture characteristics

Pseudo-dynamic output parameters

Characterization of Spatial Slip Distribution
The variability of slip in earthquakes maps directly into high frequency strong groundmotion [Bernard and Madariaga, 1983; Spudich and Frazer, 1984]. At present, it is notpossible to anticipate earthquake slip distributions deterministically, but we can anticipate thecharacteristics of slip slip variation in earthquakes of a given size if we have a valid stochastic
description of the source. With this statistical description of the earthquake source in hand, we
can generate many realizations of the slip distribution in a scenario earthquake of interest (withaccompanying strong ground motions) without the need to know the slip distribution
deterministically.
Here we use a spatial random-field model for earthquake slip to represent the slip
distribution [Mai and Beroza, 2002], but other models characterizing spatially variable slip are
possible [Andrews, 1980; Herrero and Bernard, 1994; Zeng et al., 1994; Sommerville et al.,1999]. A spatial random field is characterized either in space by its autocorrelation function,
C(r), or in the spectral domain by its power spectral density, P(k), where k is the wavenumber.
Analyzing a set 44 finite-source rupture models, Mai and Beroza [2002] found that a vonKarman auto-correlation function best represents the power spectral decay of complex
earthquake slip, with correlation length increasing with increasing magnitude. The power
spectral density P(k) of the von Karman auto-correlation function is given by

P(k) =
2

2

2

2

4 !H
ax az
H ,
K H (0) (1+ k 2 )(1)with k = ax k x + az k z , H being the Hurst exponent, and KH the modified Bessel function of
the first kind (order H). The characteristic scale lengths are given by the correlation lengths, ax
and az in the along-strike and down-dip directions. The Hurst exponent, H, in the expression for
the von Karman distribution describes the spectral decay at high wavenumbers. The spectral

6decay of 44 source models studied by Mai and Beroza [2002] revealed that H = [0.8 –1.0], andthat the correlation length scale with magnitude as:ax ≈ 2.0 + ⅓ Lefflog(ax) ≈ −2.5 + ½ Mw
(2)az ≈ 1.0 + ⅓ Wefflog(az) ≈ −1.5 + ⅓ Mw

where Leff, Weff are the effective source dimensions as defined in Mai and Beroza [2000].
We generate heterogeneous slip distributions for the selected target earthquake magnitudes
using equations (1) and (2). Source dimensions for these earthquakes (6.4 < M < 7.2) arecalculated following Mai and Beroza [2000]. Having calculated the power-spectral density P(k),the two-dimensional slip function is obtained by assuming a uniform-random phase andsubsequent two-dimensional Fourier transformation under the requirement of Hermitiansymmetry to ensure a purely real valued slip-function.

Calculating Static Stress Drop from Slip Distribution
The complete stress-time history can only be calculated by knowing the entire slip-time
evolution of the rupture [e.g., Bouchon, 1997]; however, the goal in this study is to define adynamically consistent slip-time evolution based on the final slip only. For this purpose we use
the static stress drop, Δτ, associated with slip occurring on the rupture plane, as a fundamental
parameter for characterizing the temporal slip evolution. We relate slip and stress to one anotherby a convolutional integral [Andrews, 1980] expressed as a multiplication in the wavenumber
domain, as:Δτ(k) = − K(k) · D(k)(3)

where τ(k) and D(k) denote the Fourier transforms of the two-dimensional stress drop Δτ(x,z) andslip-function D(x,z), respectively. K(k) represents the static stiffness function that for crustal
rocks can be approximated as K(k) = −½ µ k. Hence, having simulated a two-dimensional slip
distribution D(x,z), equation (3) can be used to compute the associated static stress drop Δτ(x,z)
by means of Fourier transformation of D(x,z) and subsequent inverse Fourier transformation ofτ(k).

Characterization of Fracture Energy for Temporal Slip Evolution
The characterization of fracture energy Gc is the critical step in our PD procedure fordetermining the temporal evolution of slip. Our PD approach is built on the following
relationship derived by Andrews [1976] for a simple homogeneous anti-plane crack:
1 - v2/β 2 = π2 · (Rc / 2)2,(4)

where Rc is the dimensionless parameter:
Rc = µ ·Gc / (Δτ 2·Lh),(5)

v is the rupture speed, β is the shear wave speed, µ is the shear modulus, Δτ is the stress drop, and
Lh is crack length. Guatteri et al. [2003] showed that relationship (4) applies approximately to 3-

7
D heterogeneous dynamic models suggesting that rupture velocity variability can be inferred fromthe distribution of dynamic parameters over the fault plane.
Rupture Velocity. The first step towards a temporal slip evolution characterization is todetermine a distribution of rupture velocity that is physically consistent with a given slip
realization. Given a slip distribution, Equation (4) allows us to achieve this by determining thecorresponding stress drop distribution (Equation 3), calculating the crack length Lh which we takeas the distance of each point on the fault to the hypocenter, and based on that, generating aphysically consistent distribution of fracture energy, Gc.
Fracture Energy. Different distributions of fracture energy may be consistent with a given slip
distribution and rupture propagation. Through a trial and error procedure Guatteri et al. [2003]
imposed the fracture energy distribution on their stochastic-dynamic models such that theresulting rupture velocity would be sub-shear on average everywhere over the fault plane. In
order to avoid this time consuming trial and error approach, Guatteri et al. [2002] used their setof stochastic-dynamic models to build an empirical relationship to generate a distribution of
fracture energy from a given slip model and hypocenter location. After some data exploration
they developed an empirical relationship based on the normal linear model where the conditional
expectation of Gc is:E(Gc|β, Δτ, Lh) = β0 + β1 Δτ Lh1/2,(6)where β is the vector containing β0 and β1, the regression coefficients that they determinedthrough a least squares regression procedure.Figure 3 shows the fracture energy plotted as afunction of the respective predictor Δτ Lh1/2 (that is proportional to the stress intensity factor) forall points over the fault plane for all the different dynamic models developed by Guatteri et al.
[2002].
Relationship (6) can be used either to assign a fracture energy distribution to compute
spontaneous rupture models and to apply Equation (5) to calculate a distribution of rupture
velocity in our PD procedure. Before using it for the latter application, we validated and modified
it as described below.

Figure 3: Fracture energy values as afunction of stress drop and crack
distance. Values relative to different
dynamic models described in Guatteri etal. [2003] are shown in different colors.

We first applied Equation (6) to derive a distribution of fracture energy to compute a new setof stochastic-dynamic models. However, unlike Guatteri et al. [2003], in this study we allow the
dynamic models to have areas of negative stress drop (Figure 1). This difference implies that on

8average we need lower values of Gc to maintain a sub-shear rupture velocity over the fault plane,
because the negative stress drop areas increase the resistance of the fault to rupture. As a result,
we find that lower values of regression coefficients β0 and β1 are needed than those foundpreviously by Guatteri et al. [2002].
Furthermore, we found that two distinct sets of coefficients are needed, corresponding torupture models with M <= 6.5 and for M > 6.5, respectively. The two sets are β0 = 0.18 and β1 =0.0015 for M <= 6.5, and β0 = 2.7 and β1 = 0.0021 for M > 6.5. The difference β0 stems from thedifference in average stress drop values, while the almost identical slope β1 indicates that thedependency of the stress intensity factor is preserved. We interpret this result as a consequenceof lower stress drop values in our rupture models with M <= 6.5 than those with a larger
magnitude, with corresponding different values of available elastostatic energy Eel (Figure 1,
Table 1). This would imply that the coefficients of the fracture energy distribution in Equation (6)vary as a function of earthquake magnitude.
The possibility that the scaling coefficients for fracture energy Gc are different forearthquakes below and above M = 6.5 has interesting implications as this is the magnitude range
where standard scaling relations may break down [Shimazaki, 1986; Hanks and Bakun, 2003].
Once rupture width approaches the thickness of the seimogenic zone, and rupture length L
becomes larger than about twice rupture width W, ruptures grow in along-strike direction only,
and the above defined crack length Lh becomes progressively larger for larger/longer earthquakes,
requiring higher values of fracture energy Gc
Equation (6) is derived for points on the fault having a positive stress drop. For areas ofnegative stress drop, we assign a constant value of Gc calculated as the area shown in Figure (1c)by assuming minimum allowed strength-excess and Dc values (Table 1). Note that the use of
Equation (6) provides a very simplistic parameterization of fracture energy given a slip
distribution (and therefore a static stress drop distribution) and hypocenter. The resulting fractureenergy distribution represents only one possible realization out of the space of physicallyplausible fracture energy distributions consistent with a given slip model. However, Equation (6)contributes to the intent of this paper to provide a simple tool to quickly generate a physicallyconsistent rupture model.
Having assumed in our dynamic models distributions of Dc and strength excess that are fairly
uniform and within assigned limits, the strong dependency of Gc on stress drop in our modelingappears to be an artifact of our assumptions. As pointed out by Favreau and Archuleta [2003],
however, it is plausible that larger values of stress drop are necessary to break stronger barriers ofenergy, or else the rupture would stop. The dependency found between Gc and the stress intensity
factor may be physically interpreted as corresponding to the increase of energy lost due to theoccurrence of off-fault microcracking with rupture propagation distance [Andrews, 1976; Peck etal., 1985].

Slip-Velocity Function Parameterization
A realistic parameterization of the slip-velocity function (SVF) is a critical component ofearthquake rupture modeling for strong motion prediction. As for the rupture velocity, thevariability of the SVF is poorly constrained from waveform inversion procedures. In dynamic
modeling the shape of the SVF at different fault points depends on several factors, such as local
stress parameters and distance from the hypocenter [Day, 1982; Nakamura and Miyatake, 2000;
Guatteri and Spudich, 2000; Guatteri et al., 2003]. Typical SVFs of dynamic models governedby a slip-weakening friction law are shown in Figure 4. Despite the considerable complexity seenin these slip-velocity functions, their time-dependence shows an approximate t-1/2–decay,
consistent with the quasi-dynamic Kostrov-type slip-velocity function proposed by Archuleta andHartzell [1981]. Note that some points on the fault are characterized by a long tail of low slip-

9velocity. From inversion of strong-ground motion, it is not clear whether such low-amplitude tails
would be resolved.

Figure 4: Slip-velocity functions
for identical points on the fault(each column) at a depth of 9.5 km
for four spontaneous-dynamic
rupture models [Guatteri et al,2003].

In this study, we define the rise time as the time from 5% to 95% of the total slip, i.e. being
representative of the total slip duration. A comparison between the rise time inferred for realearthquakes from kinematic inversion procedure [Wald et al. 1994; Beroza and Spudich, 1988;
Cotton and Campillo, 1995] and the rise time defined for dynamic SVF’s may not be very
meaningful. In order to characterize the shape of the SVF (impulsive vs. smooth), we define anadditional time parameter, the ‘time pulse’ Tp, as the time from 5% to 50% of the total slip. Forvery spiky SVF, Tp will be much shorter than for a broad and smooth SVF, while the respectiverise time may be of comparable size. We believe that Tp is a more relevant time parameter forstrong-ground motion prediction than the rise time, particularly at high frequencies
In this study we propose a very simple approximation of the dynamic SVF using a small
number of parameters in order both to provide a feasible parameterization within the PD
procedure and to capture the relevant characteristics of the dynamic SVF for strong-groundmotion prediction. Nakamura and Miyatake [2000] proposed a SVF parameterization that closely
resembles the SVF obtained in dynamic modeling; their characterization is perhaps superior to
our proposed approximation, but it is also rather complex to implement within our approach. Day[1982] also proposed very insightful relationships between source parameters describing the SVF
in dynamic models. We make use of both the Day [1982] and the Nakamura and Miyatake [2000]
results in formulating our parameterization. Table 2 lists the parameters defining the SVF at agiven point on the fault: the maximum slip-velocity value Vmax, the time pulse Tp, and the risetime τr.
Figure 5 shows our proposed parameterization of the SVF as a simple approximation of the
SVF typical of dynamic models governed by a slip-weakening friction law (Figure 4). The SVF is
composed of two overlapping triangles, T1 and T2 with a total base equal to the rise time, τr. T1
is an isosceles triangle having an area A equal to half of the local total slip S, height equal to Vmax
and base equal to the pulse width Tp. T1 is the portion of the SVF that contributes the most to theseismic radiation. T2 is a rectangular triangle with height equal to V2 = cVmax. The area of the
non-overlapping part of T2 is equal to A. In this study we set c = ½, but other values may be10chosen. For simplicity, let us assume that slip starts at t = t1 = 0. By solving the following
system:$* + 2Vmax 't + 2Vmax!cVmax = (T p %&)#!t = t2"(7)

t 2 = T p (1 ! c 2 ),(8)

we find t2:

and for V2 = ½ Vmax it follows that t2 = ¾ Tp.
The following system of equations applies to our SVF parameterization (Figure 5):$ 1 Vmax T p = A! 2#!" A = 1 2 V2 (& r % t 2 )% 1 2 V2 T p % t 2

(

)

(9)

Vmax
cVmax

T1

T2
t2

Tp

Figure 5. Slip-velocity function
parameterization in our PD approach,characterized by two overlapping triangles,
an isosceles triangle T1 and a rectangulartriangle T2.

!rVmax. The value of peak slip-velocity is mainly controlled by the local value of stress drop, of Dc,
and distance from the hypocenter [Day, 1982; Andrews, 1985; Ohnaka and Yamashita, 1991;
Nakamura and Miyatake, 2000; Guatteri and Spudich, 2000]. Because for a given stress drop
distribution, different choices of Dc result in different peak-slip velocities (a small Dc results in alarge Vmax, and a large Dc in a low Vmax) [Guatteri and Spudich, 2000], we cannot provide a
unique parameterization for Vmax for a given slip distribution. However, the user has the option to
set a value of Dc within a given range and tune the SVF shape in a manner consistent with thefrequency band of interest in ground motion simulation.
The analysis of our set of dynamic models suggests that the following relationship, modified
from Day [1982], for Vmax provides an adequate description of the Vmax distribution observed indynamic models:

Vmax = 0.5Vmax Re f W" p v (!µ ),(10)

where W is the width of the fault, τp = max(Δτ, w(τs − τd)) and VmaxRef = 0.9 Dcfc, where Dc is achosen value of slip-weakening distance, τs − τd = 2Gc/Dc, and fc is a frequency parameter
defined in Ohnaka and Yamashita [1989]. The equation for VmaxRef has been adapted from

11
Ohnaka and Yamashita [1989]. Our definition of τp is based on the idea that, if the value ofstrength excess is very high compared to the stress drop, then it has a large effect on Vmax. In ourparameterization, we found that w=0.6 for areas with Δτ >0 and w=1 for areas with Δτ <= 0gives a satisfactory fit with the Vmax distribution of the dynamic models.
Pulse width Tp. As Figure 4 shows, short values of Tp correspond to large values of Vmax. Figure
6 displays the relationship between Tp and the value of total slip S divided by Vmax for the points
on the fault with a slip and slip-velocity larger than 50% of the maximum slip and maximum slipvelocity over the fault, respectively. A linear relationship provides a simple parameterization of a
Tp given Vmax and the total slip:

T p = !1 S

,
(11)
Vmax
where β1 =0.84. Let Vave be the average slip-velocity during Tp, then 0.45 S = Vave Tp, implyingthat Vave Tp ≈ 1/2·Vmax Tp that justifies our parameterization 0.5·S = 0.5·Vmax Tp. Note that thisrelationship is derived for areas of large slip and slip-velocity that are the areas of the fault withlarge seismic energy radiation, and hence contribute the most to strong-ground motions.

Figure 6. Tp values of the set of 9 dynamic
models plotted as a function of thecorresponding value of slip/ V max. This plotshows the relationship for the points on thefault having slip and Vmax values larger than
50% of the maximum slip and slip-velocity,respectively, over the fault.

Slip Rise Time. As for the characterization of the fracture energy distribution, Guatteri et al.
[2003] built an empirical relationship for slip rise time based on the dynamic models, consistent
with our stochastic-dynamic models. Assuming a slip-weakening friction model [Day, 1982,
Andrews, 1985], the duration of slip at a given point on the fault is controlled mainly by the total
fault rupture duration. Figure 7 shows the rise time values τr plotted as a function of thedifference between the total effective fault rupture duration, Totrup, and the corresponding rupture
time value, Trup. The total effective fault rupture duration is the maximum rupture time value
among the fault boundary locations aligned with the hypocenter. Figure 7 suggests a normal
linear model for an empirical relationship for rise timeE(τr|β, Totrup, Trup) = β0 + β1( Totrup- Trup).

(12)

Although we have derived an empirical relationship for τr (Equation 12), we require that τr be
consistent with Equation (9). Tp is one of the most important time parameters in terms of its

12
influence on the calculated ground motion. Therefore, we follow a 2-step approach to assign all
the parameters for SVF such that the consistency with the given slip distribution is satisfied. First,
we calculate an initial guess for τr from the empirical relationship for the purpose of imposing an
upper bound to Tp. Then, from Equation (10), we can solve for τr as a function of Tp as follows:! r = Tp +

S
.
cVmax

(13)

Figure 8 shows a comparison of PD SVF’s derived by applying the parameterization describedabove with the corresponding dynamic SVF’s.

Figure 7. Slip rise time as a function of total fault
rupture duration and rupture time. Values
corresponding to different dynamic models areshown with different colors.

Figure 8. Comparison between the dynamic SVF with the PD SVF parameterization proposed in this paper for two
selected points on the fault for two different models. The rupture times of the dynamic model has been aligned to thatof the PD model for better comparison.

13

Summary of PD-approach:The basic steps involved in the development of a pseudo-dynamic source realization are thefollowing
1. Define the target earthquake magnitude and event type.
2. Assign or determine the fault dimensions [Wells and Coppersmith, 1994; Somerville etal, 1999; Mai and Beroza, 2000].
3. Set the hypocenter location.
4. Generate a stochastic slip distribution [Andrews, 1980; Somerville et al, 1999; Mai andBeroza, 2002].
5. Compute the corresponding stress drop distribution [Andrews, 1980].
6. Generate a distribution of fracture energy (Gc) using Equation (6).
7. Calculate rupture velocity and rupture time distributions using Equation (5).
8. Calculate Vmax distribution from Equation (10).
9. Apply 2-step procedure to derive the distribution of Tp and τr using Equations (11-13).
Figure 9 shows an example of a pseudo-dynamic source model for a M=6.5 strike-slip
earthquake. Figure 10 shows a comparison of a fully dynamic source model and its corresponding
pseudo-dynamic source characterization. The common constraints to the two source
characterizations are the fault geometry, the hypocenter location, the slip distribution and therequirement of sub-shear rupture propagation.

Figure 9. Example representing the pseudo-dynamic
procedure. The starting point is a slip realization generated
as a spatial random field [Mai and Beroza, 2002]. The
corresponding stress drop distribution is computed using themethod of Andrews [1980].

Figure 10. Comparison between a fully
dynamic rupture model (left) and acorresponding pseudo-dynamic source model(right). Notice that the main characteristics ofthe dynamic rupture are represented in the PD
model, such as the variation in the rupture
velocity, the areas with large and low peak slipvelocity, and the areas with short Tp.

14Discussion: The energy budget of earthquake rupture
The pseudo-dynamic source characterization developed here is not intended to be a perfect fit
to the corresponding spontaneous dynamic rupture models, which are themselves only models ofwhat actually happens in nature during earthquake rupture. The features of the earthquake source
most important for realistic strong-ground motion simulations, however, are well described by
our PD parameterization. The main correlation between source parameters can be described inplain English as follows: the distribution of stress drop is the main factor affecting the distributionof the other dynamic rupture parameters. The rupture velocity correlates with stress drop, Vmax
correlates with stress drop, and Tp and τr are inversely proportional to stress drop. Finally, theposition of the hypocenter affects the distribution of parameters, and certain hypocenter positions
are not plausible as those would not lead to spontaneous rupture propagation. To further
investigate this last statement, we analyzed the energy budget during earthquake rupture.

Plausible Models Based on Energy Budget and Hypocenter Location
It is a common and frustrating experience of many dynamic modelers to initiate spontaneous
rupture calculations that subsequently abort before rupturing to the desired earthquake size
[Nielsen and Olsen, 2000; Oglesby, 2002; R. Archuleta, pers. comm.; R. Harris, pers. comm.].
While the investigation of the underlying physical process distinguishing small and large
earthquakes is well beyond the scope of this paper, we focus on identifying those target slip
distributions that are not consistent with a given hypocenter location based on a simple energy
budget calculation. Our approach helped us both to speed up the computation of successful
spontaneous rupture models and to select physically plausible (realizable) PD models.
In describing the energy budget during earthquake rupture, we follow Favreau and Archuleta[2002]. Each point on the fault provides a seismic energy Ers = Eel – Efr – Erx, where the various
terms are defined in Table 1. For simplicity, in this study we neglect the relaxation work Erx spentat the arresting the rupture (Figure 1). As noted by Favreau and Archuleta [2002], the fault canbe characterized by locally negative seismic energy density values, but its integral on the fault
must be positive. Based on this physical requirement, we consider that, as a fundamental
condition for its growth, the rupture must be propagating such that the integral of Ers on therupture area is remains always positive. In other words, while there can be local sinks of energy
over the fault, an earthquake process cannot be described overall as a sink of energy. For a given
earthquake magnitude, we might expect a more stringent constraint such as a minimum value ofradiated seismic energy, but we leave this as an open question for future studies.
It is important to note that the placement of the hypocenter exerts a strong influence on theenergy budget as it affects both the specific fracture energy parameterization, as well as thedetermination of the areas of the fault over which we incrementally integrate Ers. The distributionof elastostatic energy Eel can be calculated from the target slip model and the corresponding stressdrop. The fracture energy distribution is assigned through Equation (6), and finally the seismicenergy density distribution is calculated as shown in Table 1. Although our calculated energy
budget strongly depends on our fracture energy parameterization, it provides a consistent
approach with our dynamic source modeling and PD procedure.
The integrated seismic energy over the ruptured fault area should remain positive in order toallow rupture growth. In other words, if the integrated seismic energy becomes negative, there is
no more energy available for the rupture to grow, and hence the rupture would stop. This is theidea proposed originally by Husseini [1977] and is consistent with the idea that a small
earthquake is like a large earthquake that ran out of energy.
Figure 11 shows the energy budget calculation for a given slip distribution (Figure 11a; thehypocenter location is shown by the star). Figure 11e shows the seismic energy integrated overconcentric areas around the hypocenter in order to mimic simplified rupture propagation. Noticethat around the hypocenter there is an area that is a large sink of energy. In Figure 11* the15corresponding spontaneous rupture model with identical stress drop and Gc distributions is
consistent with the lack of available energy necessary to maintain a rupture growth beyond thenucleation area where the rupture is forced. Figure 12 shows the energy budget calculation forthe same slip distribution but a different hypocenter. This model results in a positive integrated
seismic energy over the entire fault and is consistent with the rupture propagating over the entire
fault (Figure 12*). The same result may be achieved by lowering the Gc, and this approach could
be used to test the space of Gc and hypocenters consistent with a given slip distribution, especiallybefore running spontaneous rupture models.

Figure 11. Simulated slip distribution for a M=7
strike-slip earthquake (top left) and thecorresponding energy calculations for the PDmodel. Note of Ers (bottom left) indicatedincomplete rupture propagation.

Figure 11*. Dynamic model for the slip
distribution shown in Figure 11. Note how the
rupture terminated at an early stage, in accordance
to what Figure 11 indicates.

Figure 12. Same slip distribution
Figureas shown
12*. Dynamic
inmodel for the rupture
Figure 11, but different hypocenterparametrization
location; theshown in Figure 12. Note that the
PD-energy calculations indicatedynamicthat therupture
ruptureindeed propagated over almost
would propagate farther in this case.
the entire fault, in agreement to Figure 12, but
contrasting to Figure 11.

16Conclusions and Caveats
Dynamic rupture modeling has the advantage with respect to kinematic modeling ofproviding a physically self-consistent earthquake source characterization. Guatteri et al. [2003]
showed that it leads to a realistic representation of ground motion time histories and to a realistic
prediction of ground motion intensity from future earthquakes. In our work we have shown that a
pseudo-dynamic source characterization has the potential to improve the source model design(physically based) in probable earthquake scenarios. Pseudo-dynamic modeling is a feasible
approach for the simulation of suites of ground motion time histories, merging naturally with theprobabilistic approach usually taken for seismic hazard analysis. However, the pseudo-dynamic
procedure outline above has focused on strike-slip earthquakes in the magnitude range from 6.4
to 7.2. Future work will be devoted to investigate and develop a PD-model that also includes dipslip earthquakes, and is potentially applicable over a wider magnitude range. At the time ofwriting this report, we have only limited experience in applying the PD-approach to largemagnitude earthquakes [Mai et al, 2001] that already suggest that at least size-dependent
relationships for the fracture energy Gc should apply. The application of our PD-model toearthquake scenarios that are clearly out of the range for which the current characterization has
been developed warrants extra care and validation of the results.

References
Abrahamson, N. A. and W. J. Silva (1997). Empirical response spectral attenuation relations for shallow crustal
earthquakes, Seism. Res. Lett. 68, 94-128.
Aki, K., Scaling law of seismic spectrum, J. Geophys. Res., 72, 1217-1231, 1967.
Andrews, D.J. (1976). Rupture propagation with finite stress in antiplane strain, J. Geophys. Res., 81, 3575-3582.
Andrews, D.J. (1980). A stochastic fault model, I. Static case, J. Geophys. Res., 85, 3,867-3,877.
Andrews, D.J. (1985). Dynamic plane-strain shear rupture with a slip-weakening friction law calculated by a boundary
integral method, Bull. Seismol. Soc. Am., 75, 1-21.
Andrews, D.J. (1999). Can we infer slip-weakening displacement from seismic observations?, Eos Transactions,
American Geophysical Union, vol. 80
Arculeta, R.J., and S. Hartzell (1981). Effect of fault finiteness on near-source ground motion, Bull. Seismol. Soc. Am.,71, 939-957.
Bernard, P., and R. Madariaga, A new asympototic method for the modeling of near-field accelerograms, Bull.
Seismol. Soc. Am., 74, 539-557, 1984.
Beroza, G. C., and P. Spudich (1988). Linearized inversion for fault rupture behavior: application to the 1984, Morgan
Hill, California, earthquake, J. Geophys. Res., 93, 6275-6296.
Boatwright, J. and H. Quin (1986). The seismic radiation from a 3-D dynamic model of a complex rupture process, Part
1. Confined rupture, in Earthquake Source Mechanics, S. Das, J. Boatwright, and C. Scholz ed., AGU,
Washington, D.C.
Boore, D. M., and W. B. Joyner (1997). Site amplifications for generic rock sites. Bull. Seism. Soc. Am. 87, 327-341.
Bouchon M (1997). The state of stress on some faults of the San Andreas system as inferred from near-field strong
motion data, J. Geophys. Res., 102 (B6), 11731-11744.
Cotton, F., and M. Campillo (1995). Frequency domain inversion of strong motions: application to the 1992 Landers
earthquake, J. Geophys. Res., 100, 3961-3975.
Das. S., and B.V. Kostrov (1987). On the numerical boundary integral equation method for three-dimensional dynamic
shear crack problems, J. Appl. Mech., 54, 99-104.
Day, S.M. (1982). Three-dimensional simulation of spontaneous rupture: The effect of nonuniform prestress, Bull.
Seismol. Soc. Am., 72, 1,881-1,902.
Gelman, A., J.B. Carlin, H.S. Stern, and D. B. Rubin (1995). Bayesian data analysis. Chapman & Hall, London.
Favreau P, and R.J. Archuleta (2003). Direct seismic energy modeling and application to the 1979 Imperial Valley
earthquake, Geophys. Res. Lett., 30 (5), 1198-2002.
Guatteri, M., and P. Spudich (2000). What can strong motion data tell us about slip-weakening fault-friction laws?,
Bull. Seismol. Soc. Am., 90, 98-116.

17
Guatteri, M., P. M. Mai, and G. C. Beroza (2002). Dynamic and pseudo-dynamic source characterization for strong
ground motion prediction, Eos Trans. AGU, 83(47), Fall Meet. Suppl., Abstract S12B-1210.
Guatteri, M., P. M. Mai, G. C. Beroza and J. Boatwright (2003). Strong ground motion prediction from stochasticdynamic source models, Bull. Seismol. Soc. Am., 93, 301-313.
Hanks T.C., and W.H.Bakun (2003). A bilinear source-scaling model for M-log A observations of continental
earthquakes, Bull. Seismol. Soc. Am., 92 (5), 1841-1846.
Haskell, N., Total energy and energy spectral density of elastic wave radiation from propagating faults, II, a statisticalsource model, Bull. Seismol. Soc. Am., 56, 125-140, 1966.
Hisada, Y. (2000). A theoretical omega-square model considering the spatial variation in slip and rupture velocity, Bull.
Seism. Soc. Am. 90, 387-400.
Hisada, Y. (2001). A theoretical omega-square model considering the spatial variation in slip and rupture velocity Part
2: Case for a two-dimensional source model, Bull. Seism. Soc. Am. 91, 651-666.
Husseini, M.I. (1977). Energy-balance for motion along a fault, Geophys. J. Roy. Astr. Soc., 49 (3), 699-714.
Husseini, M.I., D.B. Jovanovich, and M.J. Randall (1975). Fracture energy of earthquakes, Geophys. J. Roy. Astr. Soc.,
43 (2), 367-385.
Ide, S. (2002). Estimation of radiated energy of finite-source earthquake models, Bull. Seism. Soc. Am. 92, 2,994-3,005.
Mai, P. M., and G. C. Beroza (2000). Source scaling properties from finite-fault-rupture models, Bull. Seism. Soc. Am.
90, 604-615.
Mai, P.M., M. Guatteri, and G. C. Beroza (2001). A stochastic-dynamic earthquake source model for strong-motion
prediction: earthquake scenarios on the Hayward-Rodgers-Creek fault system, Eos Transactions, American
Geophysical Union, 82, 869.
Mai, P. M., M. Guatteri, G. C. Beroza, and J. Boatwright (2001). Toward a more physical basis for strong-motion
simulation, Seis. Res. Lett, 72, 273.
Mai, P. M., and G. C. Beroza (2002). A spatial random-field model for complex earthquake slip, J. Geophys. Res.,
107(B11), 2308, doi:10.1029/2001JB000588.
Nakamura, H., and T. Miyatake (2000). An approximate expression of slip velocity time functions for simulation ofnear-field strong ground motion, Zishin (Journal of the Seismological Society of Japan), 53, 1-9, (in Japanese)
O’Connel, D.R., and J.P. Ake (1995). Ground motion analysis for Hoover Dam, Seismotectonic Report, 94-1, Bureau
of Reclamation, Denver, Colorado.
Ohnaka, M., and T.Yamashita (1989). A cohesive zone model for dynamic shear faulting based on experimentallyinferred constitutive relation and strong motion source parameters, J. Geophys. Res., 94 (B4), 4089-4104.
Peck, L., R.C. Nolen-Hoeksema, C. C. Barton, and R. B. Gordon (1985). Measurement of the resistance of imperfectly
elastic rock to the propagation of tensile cracks, J. Geophys. Res. 90, 7827-7836.
Quin, H.R., and S. Das (1989). A hybrid boundary integral equation method for the computation of source timefunctions for 3D rupture propagation, Geophys. J. Roy. Astron. Soc., 96, 163-177.
Shimazaki, K. (1986). Small and large earthquakes: The effects of the thickness of the seismogenic layer and the freesurface, Earthquake Source Mechanics, AGU Geophys. Mono. 37, eds S. Das, J. Boatwright, C. Scholz,
Washington D.C., American Geophysical Union, 209-216.
Somerville, P. G., N. F. Smith, R. W. Graves, and N. A. Abrahamson, (1997). Modification of empirical strong groundmotion attenuation relations to include the amplitude and duration effects of rupture directivity, Seism. Res.
Lett. 68, 199-222.
Spudich, P., and L.N. Frazer, Use of ray theory to calculate high-frequency radiation from earthquake sources having
spatially variable rupture velocity and stress drop, Bull. Seismol. Soc. Am., 74, 2061-2082, 1984.
Wald, D.J., T.H. Heaton, and K.W. Hudnut (1996). The slip history of the 1994 Northridge, California earthquake
determined from strong motion, teleseismic, GPS, and leveling data, Bull. Seismol. Soc. Am., 86, 49-70.
Wells, D.J., and K.J. Coppersmith (1994). New empirical relationships among magnitude, rupture length, rupture
width, rupture area, and surface displacement, Bull. Seismol. Soc. Am., 84, 974-1002.
Yamashita, T., and M. Ohnaka (1991). Nucleation process of unstable rupture in the brittle regime – a theoretical
approach based on experimentally inferred relations, J. Geophys. Res., 96 (B5), 8351-8367.