Research Paper

Earthquake Spectra
Simulation of near-fault 1–24
Ó The Author(s) 2021
ground motions for Article reuse guidelines:
[Link]/journals-permissions

randomized directivity DOI: 10.1177/87552930211030941

[Link]/home/eqs

parameters

Yara Daoud1, Mayssa Dabaghi1, and

Armen Der Kiureghian2,3

Abstract
The Dabaghi and Der Kiureghian stochastic near-fault ground motion model requires
information about the source, site, and source-to-site geometry, including directivity
parameters. Directivity parameters entail often unavailable knowledge of the rupture
geometry and hypocenter location. This article presents methods to randomize the
directivity parameters required to simulate near-fault ground motions. A first proce-
dure is proposed where only the contributing fault, earthquake magnitude, and site
location are known. Possible rupture directivity conditions are accounted for by ran-
domizing the rupture geometry and hypocenter location. For this purpose, new pre-
dictive models of the rupture geometry parameters are developed for shallow
crustal earthquakes with magnitudes between 5.2 and 7.9. To allow validation of syn-
thetic motions with NGA-West2 models, a second procedure randomizes the rup-
ture geometry and both hypocenter and site locations. Results show a general
agreement between the two methods.

Keywords
Near-fault ground motions, rupture geometry models, stochastic ground motion
model, directivity parameters, ground motion simulation, hypocenter location
Date received: 1 January 2021; accepted: 17 June 2021

1
American University of Beirut, Beirut, Lebanon
2
American University of Armenia, Yerevan, Armenia
3
UC Berkeley, Berkeley, CA, USA

Corresponding author:
Mayssa Dabaghi, American University of Beirut, P.O. Box 11-0236, Riad El-Solh, Beirut 1107 2020, Lebanon.
Email: md81@[Link]
2 Earthquake Spectra 00(0)

Introduction
Near-fault ground motions tend to be more complex and more variable than far-field
ground motions and may possess distinct characteristics, such as the rupture directivity,
hanging-wall (HW) effects, and the fling step, which may have profound influences on
structural response. The rupture-directivity effect depends on the direction, length, and
velocity of rupture propagation relative to a considered site (Somerville et al., 1997). When
the rupture front propagates toward the site at a velocity nearly equal to the shear-wave
velocity of the ground, all the seismic energy radiated from the fault rupture tends to arrive
at the site in a single, large-amplitude, short-duration pulse—the so-called forward direc-
tivity pulse—in the fault-normal direction. If the rupture propagates away from the site,
the site is in the backward directivity region and typically records a small amplitude, long-
duration ground motion in the fault-normal direction.
The assessment of seismic risk and performance of structures located near active faults
should account for the distinct near-fault characteristics of potential earthquakes. Given
the scarcity of recorded near-fault ground motions from larger magnitude events and exhi-
biting near-fault effects, synthetic ground motions can be useful in seismic assessment
studies in addition to or in place of recorded motions. Obviously, these synthetic ground
motions must capture the important characteristics of near-fault ground motions and
properly represent their natural variability for given earthquake source and site
characteristics.
Several ground motion simulation techniques have been proposed. Deterministic
physics-based (or source-based) techniques explicitly model the rupture process at the seis-
mic source, propagate the resulting seismic waves, and determine the site response
(Hartzell et al., 2005). Using a three-dimensional (3D) velocity model, these methods cap-
ture the local wave-propagation and basin effects. However, the generated synthetic
motions typically lack high-frequency components beyond approximately 1 Hz. To
address this shortcoming, high-frequency stochastic components are added to the physics-
based synthetics (Graves and Pitarka, 2016). Physics-based approaches require a precise
definition of the fault and rupture geometry and the shear-wave velocity of the medium,
which often is not available to engineers.
Dabaghi and Der Kiureghian (2014, 2017, 2018) proposed a site-based stochastic model
and simulation procedure to generate orthogonal horizontal components of synthetic
near-fault ground motions for specified earthquake source and site characteristics. This
method accounts for rupture directivity and produces pulse-like and non-pulse-like
motions in accordance with their observed proportions among recorded motions. The
model requires specification of the source (type-of-faulting F, earthquake magnitude Mw ,
depth to the top of the rupture plane ZTOR ), the shear-wave velocity of the top 30 m of soil
at the site (Vs30 ), and the source-to-site geometry (closest distance RRUP , directivity para-
meters sor d and uor f). As illustrated in Figure 1, depending on the nature of the fault, the
directivity parameter sor d describes the length s or width d of the portion of the rupture
that propagates between the hypocenter and the site in the direction of the slip, and the
directivity parameter uor f describes the angle u in a horizontal plane between the fault
rupture plane and the direction between the epicenter and the site, or the angle f in a ver-
tical plane between the fault rupture plane and the direction between the hypocenter and
the site. Parameters s and u are used for strike-slip faulting, while d and f are used for
dip-slip faulting (Somerville et al., 1997). The directivity parameters influence the prob-
ability that a ground motion is pulse-like and, if so, the amplitude and period of the pulse.
Daoud et al. 3

Figure 1. Directivity parameters (after the study by Somerville et al., 1997).

Ground-motion prediction equations (GMPEs) developed as part of the next-

generation attenuation (NGA)-West2 project (Bozorgnia et al., 2014) do not explicitly
account for the directivity effect, except those proposed by Chiou and Youngs (2014).
More recently, five directivity models were developed based on the NGA-West2 ground
motion database and numerical simulations of large earthquakes (Spudich et al., 2013,
2014). Comparisons of these models showed that their predictions have high variability
and strong dependence on individual assumptions, especially for non-strike-slip faults
(Spudich et al., 2013, 2014). These models are not considered sufficiently developed to be
incorporated into the GMPEs and, therefore, further studies have been recommended
(Spudich et al., 2013, 2014).
The near-fault simulation procedure of Dabaghi and Der Kiureghian (2018) requires as
input the parameter set (F, Mw , ZTOR , RRUP , Vs30 , sor d, uor f). To perform seismic design or
assessment at a site with known location and Vs30 , structural engineers can obtain informa-
tion about the magnitude Mw , the distance RRUP , and the source faults of earthquakes that
contribute most to the hazard from deaggregation of probabilistic seismic hazard analysis,
for example, as provided by the US Geological Survey (USGS, 2019) unified hazard tool.
However, structural engineers are not likely to have information about the parameters
sor d, uor f, or even ZTOR , which require knowledge of the rupture geometry as well as the
location of the hypocenter.
Various models have been proposed to describe the scaling of rupture geometry para-
meters (e.g. length, width, aspect ratio, area) with magnitude and vice versa (Ellsworth,
2003; Goda et al., 2016; Hanks and Bakun, 2002; Leonard, 2010; Shaw, 2009;
Thingbaijam et al., 2017; Wells and Coppersmith, 1994). These models employ various
functional forms and are developed by regression analysis using diverse historical datasets
of earthquakes from different tectonic regions. For a given Mw , they can predict different
4 Earthquake Spectra 00(0)

rupture dimensions. With the exception of the study by Goda et al. (2016), these models
predict the rupture length and width independently, without accounting for possible corre-
lation between them. Moreover, they do not account for possible differences in the scaling
of rupture geometry parameters between buried and surface ruptures.
Fewer studies are available for predicting ZTOR . Gupta (2006) used the deterministic
value ZTOR = 3 km, independent of Mw ; Campbell et al. (2009) used a median estimate of
ZTOR that depends on Mw , while Kaklamanos et al. (2011) estimated ZTOR based on hypo-
central depth, down-dip rupture width WR , and dip angle d, assuming a deterministic hypo-
center location at 60% down WR . Chiou and Youngs (2014) proposed a model that relates
the mean of ZTOR to Mw and the type of faulting, and used it in their own NGA West2
GMPE. However, their model does not account for the variability of ZTOR .
The objective of this article is to develop procedures to simulate the parameters ZTOR ,
RRUP , sor d, and uor f, when they are unknown, to generate pairs of horizontal components
of synthetic near-fault ground motions using the model and simulation procedure of
Dabaghi and Der Kiureghian (2018). Two procedures are developed as described below.
The first procedure (P1) generates the parameters (ZTOR , RRUP , sor d, uor f) randomly for a
specified site, given the type of faulting (F), the magnitude (Mw ) of the earthquake, the geo-
metry of the contributing fault (location, strike, dip, maximum dimensions), and the location
of the site. This is achieved by randomizing the fault rupture geometry, its location within
the fault plane, and the hypocenter location within the rupture plane, based on their respec-
tive probability distributions. First, the rupture geometry is simulated using distributions of
the rupture length LR , rupture width WR , and ZTOR . These distributions are developed herein
using regression analysis on a subset of earthquakes in the NGA-West2 database. Next, the
location of the rupture plane is randomized within the fault plane. The hypocenter location
along the strike and down-dip is then simulated based on the model developed by Mai et al.
(2005). Finally, the coordinates of the site and of the simulated rupture plane and hypocen-
ter are used to calculate RRUP and the directivity parameters sor d and uor f. This procedure is
useful for structural engineers, who wish to simulate sets of acceleration time series for a par-
ticular site of interest for seismic design or assessment studies.
The second procedure (P2) randomly generates sor d and uor f for given F, Mw , RRUP ,
and ZTOR and the dip angle of the fault. This procedure simulates the conditions present in
the NGA-West2 dataset and is used to compare the synthetic ground motions with the
GMPEs. This is achieved by randomizing the fault rupture geometry, the hypocenter loca-
tion within the rupture plane, and the site location. The rupture geometry and hypocenter
location are simulated as described above, while the site location is sampled from a uni-
form distribution along the locus of points at a distance RRUP from the simulated fault
rupture.
After a brief review of the Dabaghi and Der Kiureghian (2014, 2017, 2018) stochastic
near-fault ground motion model and simulation procedure, we develop predictive models
for the rupture geometry parameters LR , WR , and ZTOR in terms of the earthquake magni-
tude and type-of-faulting, and compare them to existing models. Our models are limited
to shallow crustal earthquakes in active tectonic regions having magnitudes between 5.2
and 7.9. The hypocenter location model developed by Mai et al. (2005) and adopted in this
article is then briefly reviewed. Next, procedure P1 is described and illustrated by simulat-
ing ground motions for a specified site given F and Mw and the geometry of the contribut-
ing fault. Procedure P2, which generates sor d and uor f for a random site for given F, Mw ,
Daoud et al. 5

RRUP , and ZTOR and the dip angle of the fault, is then described. Finally, near-fault ground
motions simulated using P2 are compared with the NGA-West2 GMPEs.

Stochastic near-fault ground motion model and simulation method

Dabaghi and Der Kiureghian (2014, 2017) developed a site-based parameterized stochastic
model of near-fault ground motion in two orthogonal horizontal directions. The model
accounts for the rupture directivity effect and considers both pulse-like and non-pulse-like
motions. The model employs a modulated and filtered white-noise (MFW) process with
time-varying filter parameters (Rezaeian and Der Kiureghian, 2008), and a modified ver-
sion of the Mavroeidis and Papageorgiou (2003) pulse model. The model is defined in
terms of physically meaningful parameters that represent the amplitude, duration, fre-
quency content, and pulse characteristics of near-fault ground motions. The model was
fitted to recordings taken from the Pacific Earthquake Engineering Research Center
(PEER) NGA-West2 database to obtain empirical ‘‘observations’’ of the model para-
meters. The dataset used was limited to near-fault recordings from shallow crustal strike-
slip, reverse, or reverse-oblique earthquakes in active tectonic regions with 5.5 < Mw < 8.0
and recorded at sites located at RRUP \ 31 km. Empirical predictive equations were con-
structed for the model parameters in terms of earthquake source and site characteristics,
namely F, Mw , ZTOR , RRUP , Vs30 , sor d, and uor f. The variability of the model parameters
and statistical correlations between them were estimated. Dabaghi and Der Kiureghian
(2014, 2018) extended their model to generate suites of horizontal pairs of synthetic near-
fault ground motions for specified earthquake scenarios. Figure A1 in Supplemental
Appendix A shows a flowchart of this simulation procedure, where the empirical model of
the probability of occurrence of the directivity pulse by Shahi and Baker (2014) is used to
generate pulse-like and non-pulse-like ground motions in appropriate proportions. The
resulting motions were shown to have realistic time series and to reproduce important fea-
tures of recorded near-fault ground motions, including the directivity effect and natural
variability. See the study by Dabaghi and Der Kiureghian (2014, 2018) for more details.

Rupture geometry
We wish to generate random realizations of the parameters (ZTOR , RRUP , sor d, uor f) accord-
ing to probability distributions that are consistent with the physics of the earthquake rup-
ture process and the available data. We start by developing new predictive models for
ZTOR , LR , and WR with F and Mw as the predictor variables, but consideration is also given
to the dip angle d of the fault and the seismogenic depth H (thickness of the crust), if they
are available. Differences in the scaling of rupture dimensions for buried (ZTOR .0) and
surface (ZTOR = 0) ruptures are also considered. The study is limited to crustal earthquakes
in active tectonic regions.

Database of earthquake ruptures

A subset of records from the PEER NGA-West2 database (PEER, 2015) was used to fit
the Dabaghi and Der Kiureghian (2014, 2018) model. The NGA-West2 database, which is
limited to shallow crustal earthquakes in active tectonic regions, also provides the rupture
geometry and dimensions for each earthquake. This database is used again in this study to
fit the rupture geometry models. The selected subsets are those earthquakes whose rupture
geometry and dimensions are determined by inversion of seismic waveforms. Included are
6 Earthquake Spectra 00(0)

Figure 2. Distribution of earthquakes in the considered database with respect to moment magnitude
and type of faulting.

77 earthquakes with 5:2<Mw <7:9 from strike-slip (SS), reverse (REV), and normal (N)
faulting mechanisms (see Table B1 in Supplemental Appendix B for more details). Figure 2
is a scatter plot showing the distribution of earthquakes in the considered database with
respect to moment magnitude and type-of-faulting. We note that data from normal and
normal–oblique faults as well as data from reverse and reverse–oblique faults are combined
throughout the analysis.

Depth to top of rupture, Z TOR

Data indicate that ZTOR is zero-inflated, that is, the zero value is frequently observed, as
can be seen in Figure 3. Thus, the value of ZTOR is used to classify ruptures as surface
(ZTOR = 0), that is, the rupture extends to the Earth’s surface, or buried (ZTOR .0), that is,
the rupture does not reach the surface. We note that other references define a surface rup-
ture differently. For example, Kagawa et al. (2004) define a surface rupture as having a
clear surface dislocation and significant slip, and as a buried rupture otherwise.
We employ a two-part model to represent the zero and non-zero values of ZTOR . The
first is a logistic regression model that predicts the probability that ZTOR .0. The second is
a linear regression model that predicts the mean and variability of the non-zero data.
While developing these models, it was found that the ZTOR data from normal faulting are
not statistically different from the data obtained for strike-slip or reverse faulting. This
could be due to the limited number of earthquakes with normal faults in the considered
database. For this reason and following the study by Chiou and Youngs (2014), in this
study, the data from normal faulting are combined with the data from strike-slip faulting
to fit the ZTOR models. The fitted logistic regression model is given by:
8
> 1
< , for SS and N
 ð15:2932:293M wÞ
PrðZTOR .0Þ = 1 + e 1 ð1Þ
>
: , for REV
1 + eð17:2202:293Mw Þ
Daoud et al. 7

Figure 3. Earthquake magnitude versus depth to the top of rupture in the considered database.

Figure 4. Logistic regression models for Pr(ZTOR .0) against used data.

and is plotted in Figure 4 alongside the data used for fitting.

To develop predictive equations for ZTOR when ZTOR .0, and similarly for LR and WR in
the next section, the data are first transformed to ensure normality of the regression resi-
duals. In all cases, the logarithmic transformation is used. The functional form of the pre-
dictive equation of a transformed model variable Y is initially constructed as:

Y = ðb0 + b1 F Þ + ðb2 + b3 F ÞMw + E ð2Þ

where E denotes the regression residual, assumed to be a normally distributed random vari-
able with a zero mean and standard deviation s, and b0 , ..., b3 are the model parameters.
F = 0, 1, and 2 correspond to SS, REV, and N faults, respectively. Step-wise linear
8 Earthquake Spectra 00(0)

regression is utilized to investigate the importance of the predictors, and the models are
refitted accordingly. The models are evaluated based on the resulting R2 value, residual
plots, plots of fitted versus observed values, and normality of the residuals using the
Shapiro–Wilk normality test (Shapiro and Wilk, 1965). A significance level of 5% is used
for all statistical tests unless stated otherwise.
As previously mentioned, following the study by Chiou and Youngs (2014), data from
N faults are combined with SS faults when developing the ZTOR models. Analysis of the
data for ZTOR .0 leads to the model:

9:361  1:400Mw for SS and N
E½ln ZTOR jZTOR .0 = , sln ZTOR = 0:860 ð3Þ
6:362  0:789Mw for REV

where ZTOR is measured in kilometers, and E½: and s are the expected value and standard
deviation from the regression model, respectively. Moreover, to prevent the model from
giving unreasonable predictions at lower magnitudes and in the tail, upper limits are
imposed on ZTOR . Based on observed maximum values in the database, the following
upper bounds are considered:

expðE½ln ZTOR jZTOR .0, SS and N] + 1:75sln ZTOR Þ
ZTOR < , for SS and N
 10 km ð4Þ
expðE½ln ZTOR jZTOR .0, REV + 1:75sln ZTOR Þ
ZTOR < , for REV
15 km

Figure 5 shows the predictive model for ZTOR .0 given by Equation 3 and the upper
bounds given by Equation 4, plotted against the data used for fitting. It can be seen in
Figure 4 that the probability of having a buried rupture decreases with increasing magni-
tude and is larger for reverse faults than for strike-slip or normal faults. For buried rup-
tures in Figure 5, ZTOR tends to decrease with magnitude, and reverse faults tend to have
larger ZTOR values than strike-slip and normal faults. The observed differences between
fault types could be explained by the fact that strike-slip and normal faults tend to be stee-
per than reverse faults; thus, the rupture width of a reverse fault can extend further along
the dip before rupturing the surface (or reaching the seismogenic depth). These

Figure 5. Predictive models for ZTOR .0, with upper limits, against used data.
Daoud et al. 9

Figure 6. Comparison of developed predictive models for LR with WC94 and employed data.

observations are consistent with previous ZTOR models, such as the one proposed by Chiou
and Youngs (2014). The reader is cautioned against extrapolation of these models beyond
the range of the data (5:2<Mw <7:9).

Rupture dimensions LR and W R

As mentioned before, calculating RRUP , sor d, and uor f for a given site requires knowledge
of the geometry and coordinates of the fault rupture and the location of the hypocenter.
The fault rupture is idealized as a rectangle with the aspect ratio varying from earthquake
to earthquake. In this section, predictive models are developed for LR and WR . Similar to
the study by Wells and Coppersmith (1994) and Thingbaijam et al. (2017), we differentiate
between strike-slip, reverse and normal faulting because of the differences in their rupture
mechanisms and dip angles. Possible influence on the rupture dimensions depending on
the buried (ZTOR .0) or surface (ZTOR = 0) nature of the rupture is also considered.
A procedure similar to that used for developing the linear regression model of ZTOR .0 is
used. For the rupture length, no statistical difference is observed between surface and bur-
ied ruptures, so a predictive model that does not make this distinction is developed. The
result is as follows:
8
< 5:821 + 1:402Mw for SS,
E½ln LR  = 6:106 + 1:402Mw for REV, sln LR = 0:401 ð5Þ
:
5:919 + 1:402Mw for N,

where LR is in kilometers. Note that the rate of change with magnitude is not found to be
affected by the type of faulting; only the intercept is dependent on the type of faulting.
Figure 6 shows the fitted predictive model for LR together with the employed data. The
model indicates that LR increases with Mw , and that SS ruptures tend to be longer than
REV and N ruptures. These are consistent with previous findings (Schwartz, 2018). Figure
6 also shows the widely used model for subsurface rupture length developed by Wells and
Coppersmith (1994), referred to as WC94. The comparison shows that, in general, the
developed models are similar to the WC94 models and are a good fit to the available data
for LR .
10 Earthquake Spectra 00(0)

Figure 7. Comparison of developed predictive models for WR with WC94 and employed data.

For the rupture width, a significant statistical difference is observed between the surface
and buried rupture cases. Thus, the following models are developed:
8
< 0:860 + 0:263Mw
> for SS,
E½ln WR jZTOR = 0 = 1:591 + 0:263Mw for REV, sln WR jZTOR = 0 = 0:278 ð6Þ
>
:
0:958 + 0:263Mw for N,
8
< 3:403 + 0:907Mw
> for SS,
E½ln WR jZTOR .0 = 3:211 + 0:907Mw for REV, sln WR jZTOR .0 = 0:286 ð7Þ
>
:
3:164 + 0:907Mw for N,

where WR is in kilometers. Note again that the rate of change with magnitude is not statis-
tically found to be affected by the type of faulting; only the intercept is dependent on the
type of faulting. Figure 7 shows the fitted predictive models for WR , together with the data
used to fit them. The models indicate that WR increases with Mw , and that REV ruptures
tend to be wider than SS ruptures. As for N faults, their WR scaling is similar to that of
REV faults for buried ruptures, and close to (slightly above) that of SS faults for surface
ruptures. The models also indicate that, when the rupture reaches the surface, the rate of
change of WR with Mw decreases, which is indicative of possible saturation due to the finite
width of the seismogenic layer, while the rate of change of LR with Mw is not affected (see
Equations 5 to 7). This is consistent with the observation of Dalguer et al. (2008) that sur-
face ruptures tend to have a larger aspect-ratio than buried ruptures. Again here, the larger
width of REV surface ruptures compared to SS and N surface ruptures could be explained
by the smaller dip angles of REV faults compared to the steeper SS and N faults. Figure 7
also shows the WC94 rupture width models, which distinguish between the three faulting
mechanisms but not between surface and buried ruptures.
Upper bounds on LR and WR are imposed based on seismological considerations to
avoid unrealistic values. If the length LF of the fault is known, the upper bound LR, max = LF
is imposed. Otherwise LR, max = 500km is assumed, based on the largest observed rupture
Daoud et al. 11

Table 1. Estimated correlation matrix of regression residuals and their 95% confidence intervals (in
brackets)
For ZTOR = 0 For ZTOR .0
ln LR ln WR ln ZTOR ln LR ln WR

ln LR 1 0.026 ln ZTOR 1 20.028 20.121

[20.36, 0.41] [20.30, 0.25] [20.38, 0.16]
ln WR 0.026 1 ln LR 20.028 1 0.244
[20.36, 0.41] [20.30, 0.25] [20.03, 0.49]
ln WR 20.121 0.244 1
[20.38, 0.16] [20.03, 0.49]

length (Schwartz, 2018). The upper bound WR, max depends on the width WF of the fault. If
WF is not known, WF = HZ sind is assumed, where H is the seismogenic depth, ZTOF is the
TOF

depth to the top edge of the fault, and d is the dip of the fault in degrees. If ZTOF or H are
not available, default values ZTOF = 0 and H = 25 km (Watts and Burov, 2003) are assumed.
If the dip angle is not available, mean values based on faulting mechanism may be used:
85° for SS faulting, 40° for REV faulting, and 53° for N faulting. These mean values are
obtained based on fault dip data in the considered database of shallow crustal earthquakes
and are consistent with the values reported by Thingbaijam et al. (2017). Thus, for a given
value of ZTOR , the upper bound on the rupture width is calculated as:

ðZTOR  ZTOF Þ
WR, max = WF  ð8Þ
sind
Next, correlations between the log-transformed rupture geometry parameters ( ln ZTOR ,
ln LR , ln WR ) are examined. These are estimated as the correlations between the residuals
e of the regression models (see Equation 2) and are listed in Table 1 along with their 95%
confidence intervals (in brackets; Witte and Witte, 2016). The data indicate that there is a
mild positive correlation of 0.244 (at 8% significance level) between the logarithms of the
rupture length and width only for buried ruptures (ZTOR .0), while for surface ruptures,
length and width are effectively uncorrelated. Note that when surface and buried ruptures
are combined, the correlation between the logarithms of the rupture length and width is
found to be 0.150, consistent with the correlation reported in the study by Goda et al.
(2016). Further investigation shows that the correlation between ln LR and ln WR tends to
decrease with magnitude, which explains why buried ruptures, which are usually associ-
ated with smaller magnitude events, show a positive correlation, while surface ruptures,
which are usually associated with larger magnitude events, show close to zero correlation.
Moreover, for buried ruptures, ZTOR appears to be uncorrelated with both LR and WR .
Later, for simulation purposes, we assume a correlation of 0.244 between ln LR and ln WR
for buried ruptures and zero correlation among ln LR , ln WR , and ln ZTOR for all other
cases.

Simulation of rupture plane within fault plane

Given the coordinates of the top points of the fault, the fault length LF and strike are
determined. In addition knowing the fault width WF and dip angle d, the rectangular fault
plane is defined. For any set of the simulated parameters (LR , WR , ZTOR ), the location of
12 Earthquake Spectra 00(0)

the LR 3WR rectangular rupture at depth ZTOR within the fault plane is randomized assum-
ing a uniform distribution within the fault plane.

Simulation of hypocenter location within rupture plane

The location of the hypocenter within the rupture plane has an important influence on the
near-fault ground motion since it influences the directivity effect (Somerville et al., 1997).
The location of the hypocenter is a required input for physics-based simulations (Graves
and Pitarka, 2016; Hartzell et al., 2005; Pitarka et al., 2017) and it is also required to deter-
mine the directivity parameters sor d and uor f of the Dabaghi and Der Kiureghian (2018)
model. To properly account for the variability of the generated near-fault ground motions,
the simulation procedure should account for the variability in the hypocenter location.
Some previous simulation studies, for example, Graves and Pitarka (2016) and Pitarka
et al. (2017), have neglected the variabilities in the fault rupture area and hypocenter
location.
Most previous studies assume the location of the hypocenter is uniformly distributed
along the fault strike. The down-dip location is assumed to be either deterministic or to
also follow a uniform distribution. For example, early studies in the CyberShake project
placed hypocenters every 20 km along the strike (a discrete uniform distribution) and at a
down-dip distance of 0.75 of the rupture width (Graves et al., 2011), while more recent
studies placed them every 4.5 km along the strike and down-dip (Southern California
Earthquake Center (SCEC), 2018). However, Mai et al. (2005) found that data from more
than 50 earthquakes do not support the assumption of uniform distribution of the hypo-
center location over the rupture plane. Their data indicated that the hypocenter is more
likely to be near the center of the rupture plane than near the edges. They modeled the
location of the hypocenter along the strike as a random variable following a truncated
normal distribution with a mean of 0.5LR and a standard deviation of 0.23LR, while the
down-dip location normalized by the rupture width WR was modeled using a truncated
Weibull distribution. The estimated shape and scale parameters of the fitted Weibull distri-
butions are, respectively, 0.626 and 3.921 for SS faults and 0.692 and 3.394 for REV and
N faults. We adopt these models in this article, while assuming that the along-strike and
down-dip locations of the hypocenter within the rupture plane are statistically
independent.

Procedure P1: simulation for specified site and random rupture

geometry and hypocenter location
We first consider the case where F, Mw , the fault geometry (coordinates of the top points,
fault dimensions LF and WF , depth to the top edge of the fault ZTOF , and dip angle d) and
the site location are known, that is, procedure P1. In addition, the site Vs30 is assumed to
be known. As mentioned earlier, if some of the fault geometry parameters are unavailable,
the default values WF = HZ 8
sind , ZTOF = 0 km, H = 25 km and d = 85 for SS faulting, 40° for
TOF

REV faulting, and 53° for N faulting are used. The following procedure is employed to
simulate ZTOR , RRUP , and the directivity parameters sor d and uor f:

1. Given Mw and F, compute the probability p of obtaining ZTOR .0 from Equation 1.

Randomly generate the number 0 or 1 with probabilities 1  p and p, respectively.
If number 0 is generated, ZTOR = 0, otherwise, ZTOR .0.
Daoud et al. 13

2. If ZTOR = 0, that is, the case of surface rupture, simulate ln LR and ln WR as normally
distributed uncorrelated random variables with means and standard deviations
computed from Equations 5 and 6, respectively. If ZTOR .0, that is, the case of bur-
ied rupture, simulate ln ZTOR , ln LR , and ln WR as jointly normally distributed ran-
dom variables with means and standard deviations estimated from Equations 3, 5,
and 7, respectively, and a correlation coefficient of 0.244 between ln LR and ln WR .
If any of the simulated variables exceeds its specified upper bound, repeat step 2
until all upper bounds are satisfied.
3. Given the geometry of the fault plane and the simulated rupture parameters (ZTOR ,
LR , and WR ), uniformly sample the location of the rupture plane within the fault
plane and calculate the corresponding coordinates of the edges of the rupture
plane.
4. Given F, LR , and WR , simulate the location of the hypocenter along the strike and
down-dip according to the Mai et al. (2005) model. Use the coordinates of the rup-
ture plane and d to calculate the coordinates of the hypocenter.
5. Use the coordinates of the site and the simulated coordinates of the rupture plane
and the hypocenter to determine RRUP , sor d, and uor f.

Each simulation results in a different set of ZTOR , LR , WR , and rupture and hypocenter
locations, thus in different values of RRUP , sor d, and uor f. This process is repeated to gener-
ate any desired number of rupture directivity scenarios for a specified set of input values
(F, Mw , fault geometry, and site location). In addition, given VS30 of the site, the resulting
(F, Mw , ZTOR , RRUP , Vs30 , sor d, uor f) scenarios are provided as input to the near-fault simu-
lation procedure of Dabaghi and Der Kiureghian (2018) to generate the desired number
of synthetic near-fault ground motions. Figure 8 presents simplified flowcharts of the
Dabaghi and Der Kiureghian (2018) simulation procedure (left) and of simulation proce-
dure P1 (center), and the relation between them. A more detailed flowchart of P1 can be
found in Figure A2 of Supplemental Appendix A.

Figure 8. Simplified flowcharts of the Dabaghi and Der Kiureghian (2014, 2018) near-fault (NF) ground
motion (GM) simulation procedure (left), simulation procedure P1 (center), and simulation procedure P2
(right). More detailed flowcharts of each of these simulation procedures can be found in Supplemental
Figures A1 to A3.
14 Earthquake Spectra 00(0)

Table 2. Ranges of simulated rupture geometry parameters and calculated distance and directivity
parameters
Simulated parameter LR (km) WR (km) ZTOR (km) RRUP (km) d (km) f (°)

Min 7.9 5.2 3.1 6.4 1.7 0.4

Max 18.4 20.3 9.3 10.8 17.0 16.8

Example application
The proposed simulation procedure is used to generate pairs of horizontal components of
near-fault ground motions for an earthquake event at a site selected from the CyberShake
platform (SCEC, 2018) in downtown Los Angeles (LADT) having latitude 34.05204, longi-
tude 2118.25713, and Vs30 = 390 m/s. The 2008 edition of the USGS hazard model together
with its unified hazard tool is used to determine the earthquakes and faults that contribute
most to the hazard at the selected LADT site. For example, at a spectral period of 2 s and
for a return period of 2475 years, the deaggregation results show a major contribution
from earthquakes having large magnitudes (Mw .6:50) and occurring at short distances
(\ 10 km) from the site. One of the largest contributors to the hazard is an earthquake
with Mw = 6:55 on the Upper Elysian Park reverse fault. This earthquake scenario is used
to illustrate our proposed simulation procedure. Information about the geometry of the
Upper Elysian Park fault is extracted from the CyberShake platform, which is based on
Version 2.0 of the Uniform California Earthquake Rupture Forecast, UCERF2.0 (Field
et al., 2009). The extracted information includes the coordinates of the endpoints of the
top edge of the fault, from which the fault length LF = 18:7 km is calculated as the horizon-
tal distance between the two points, that is, by idealizing the fault as a rectangular plane.
The above reference also provides WF = 20:4 km, ZTOF = 3:1 km, and d = 508.
Given information about the fault type F = 1 (REV) and geometry, and the site loca-
tion, procedure P1 is used to generate, for Mw = 6:55, 100 different realizations of ZTOR , LR ,
WR , and rupture and hypocenter locations within the fault plane, which are in turn used to
calculate the corresponding RRUP , sor d, and uor f values. Table 2 summarizes the resulting
ranges of some of these simulated variables and Figure 9 illustrates selected realizations of
the rupture geometry and hypocenter locations.
For the 100 different sets of simulated rupture realizations, the corresponding para-
meters ZTOR , RRUP , d, and f together with F = 1 (REV), Mw = 6:55, and Vs30 = 390 m/s are
used as input to the simulation procedure of Dabaghi and Der Kiureghian (2018) to
obtain 100 pairs of horizontal components of near-fault ground motions. The 5% damped
RotD50 pseudo-acceleration response spectra (Boore, 2010) of the simulated motions and
their median and median plus and minus one logarithmic standard deviations are illu-
strated in Figure 10a. These simulations, which cover a range of directivity configurations
(see Table 2), consist of 38 pulse-like and 62 non-pulse-like motions. For the former set,
the pulse period Tp ranges from 0.44 to 2.44 s. Over this period range, the pulse-like
motions tend to have larger amplitudes than non-pulse-like motions. Moreover, Figure
10a highlights one pulse-like simulation (thick dark gray line) with Tp = 2:12 s and resulting
from a forward directivity scenario (d = 11:2 km), and one non-pulse-like simulation (thick
light gray line) resulting from a backward directivity scenario (d = 3:8 km). Note that the
spectral amplitudes of the pulse-like motion are amplified near the period of the pulse.
The acceleration, velocity, and displacement time series of the two highlighted synthetic
ground motions are shown in Figure 10b.
Daoud et al. 15

Figure 9. Selected realizations of rupture geometry (gray rectangular surface) and hypocenter location
(asterisk) along strike (HypX) and down-dip (HypZ) for an Mw = 6:55 earthquake on the Upper Elysian
Park fault (dashed rectangle).

Procedure P2: simulation for random site, rupture geometry and

hypocenter location
To compare the simulation results with NGA-West2 GMPEs, we consider the case where,
in addition to F, Mw , and the Vs30 of the site, RRUP and ZTOR are also known, that is, proce-
dure P2. Unless available, the default values are used for the fault geometry parameters LF ,
WF , ZTOF , and d, similar to P1. The following procedure is developed to simulate rupture
directivity conditions (sor d, uor f):

1. Given Mw , F, and ZTOR , simulate ln LR and ln WR as normally distributed uncorre-

lated random variables with means and standard deviations estimated from
Equations 5 and 6, respectively, if ZTOR = 0, or as jointly normal random variables
with means and standard deviations estimated from Equations 5 and 7, respec-
tively, and correlation coefficient 0.244, if ZTOR .0. Follow a procedure similar to
step 2 in P1 to ensure sampled LR and WR values are within the specified upper
bounds.
2. Given F, LR , and WR , simulate the location of the hypocenter within the rupture
plane similar to step 4 in P1. Use ZTOR and d to calculate the coordinates of the
hypocenter.
3. Given the simulated rupture geometry, we determine a ‘‘racetrack’’ at a distance
RRUP from the simulated fault rupture using the method developed by Brian Chiou
(2018, personal communication). The ‘‘racetrack’’ represents the locus of all points
on the ground surface that are located at a distance RRUP from the simulated rup-
ture plane. Use the uniform distribution to sample a single site on the racetrack.
4. Use the simulated rupture geometry and hypocenter and site locations to determine
the directivity parameters sor d and uor f.
16 Earthquake Spectra 00(0)

Figure 10. (a) RotD50 spectra of 100 simulated motions (38 pulse-like and 62 non-pulse-like) at the
LADT site (Vs30 = 390 m/s) due to an earthquake with Mw = 6:55 occurring on the Upper Elysian Park
fault; median and median plus and minus one logarithmic standard deviations spectra, and examples of
one pulse-like motion (from a scenario with d = 11:2 km and with Tp = 2:12 s) and one non-pulse-like
motion (from a scenario with d = 3:8 km) are shown. (b) Acceleration, velocity, and displacement time
series of the two horizontal components of the pulse-like (dark gray) and non-pulse-like (light gray)
ground motions that are highlighted in (a).
Daoud et al. 17

Each simulation results in different LR , WR , and hypocenter and site locations and,
therefore, different values of sor d and uor f. This process is repeated to generate any desired
number of rupture directivity realizations for specified F, Mw , RRUP , ZTOR , and d. This
allows a fair comparison of the synthetic motions with the NGA-West2 ground motion
prediction equations. Figure 8 (right) presents a simplified flowchart of simulation proce-
dure P2, and its relation with the Dabaghi and Der Kiureghian (2014, 2018) simulation
procedure. A more detailed flowchart of P2 is presented in Figure A3 in Supplemental
Appendix A.

Example application
The proposed simulation procedure P2 is used to generate 600 pairs of orthogonal hori-
zontal components of near-fault ground motions for different hypothetical earthquake sce-
narios with selected F, Mw ,RRUP , ZTOR , and Vs30 values and random directivity conditions.
For each simulation, LR , WR , and hypocenter and site locations are generated following
the procedure described in the previous section. The resulting directivity variables are used
as input to the simulation procedure of Dabaghi and Der Kiureghian (2018) to obtain the
600 pairs of orthogonal horizontal components of the ground motion.
The statistics of elastic response spectra of the generated synthetic ground motions are
compared to those of the NGA-West2 GMPEs (Bozorgnia et al., 2014), and by that indir-
ectly to recorded ground motions. We use a weighted average of the five NGA-West2
GMPEs developed by Abrahamson et al. (2014), Boore et al. (2014), Campbell and
Bozorgnia (2014), Chiou and Youngs (2014), and Idriss (2014), respectively, referred to as
ASK14, BSSA14, CB14, CY14, and I14. Following the recommendation by Rezaeian
et al. (2014), all models are assigned a weight of 2/9 except I14, which is assigned a weight
of 1/9. The weighted average of the five NGA-West2 GMPEs is denoted the ‘‘NGA-
West2 model.’’
The above GMPEs have limitations in the near-fault region. First, they do not require
directivity parameters as input (except CY14) and, therefore, do not differentiate between
forward and backward directivity sites. Thus, effectively, they represent random directivity
conditions. Moreover, they were calibrated to broader ranges of magnitude and distance
than the Dabaghi and Der Kiureghian (2018) model, which is limited to records from large
magnitude earthquakes (Mw .5:5) at short distances (RRUP\31 km). Therefore, the NGA-
West2 GMPEs are not specifically calibrated to near-fault scenarios.
The synthetic ground motions are compared with the GMPEs for several hypothetical
earthquake scenarios occurring on vertical SS (F = 0) or dipping REV (F = 1) faults, having
magnitudes Mw = 6.5, 7, or 7.5, and ZTOR = 0 or 3 km. The considered sites have Vs30 =
360, 525, or 760 m/s and are located at RRUP = 5, 10, or 20 km from the fault rupture.
Simulations were conducted for a total of 34 scenarios that are reported in the study by
Daoud et al. (forthcoming). For each of the considered scenarios, statistics of the 5%
damped pseudo-acceleration response spectra of the 600 synthetic motions are compared
with those described by the NGA-West2 model. The compared statistics are the median
and median plus and minus one logarithmic standard deviation levels for the RotD50 hori-
zontal component.
The NGA-West2 GMPEs require additional input parameters not required by the
Dabaghi and Der Kiureghian model. For instance, some GMPEs need as input other dis-
tance parameters in addition to or instead of RRUP . These distance parameters include the
18 Earthquake Spectra 00(0)

Figure 11. Median and median plus and minus one logarithmic standard deviation of RotD50 spectra of
600 synthetic motions; median spectra of subsets of pulse-like (P) and non-pulse-like (NP) synthetic
motions; median and median plus and minus one logarithmic standard deviation spectra predicted by a
combination of five NGA-West2 GMPEs, and median spectra predicted by each of the five NGA-West2
GMPEs for two different earthquake scenarios.

Joyner–Boore distance Rjb (i.e. the closest distance to the surface projection of the coseis-
mic rupture), the horizontal distance Rx from the top of the rupture measured perpendicu-
lar to the fault strike, and the horizontal distance Ry0 off the end of the rupture measured
parallel to the strike. Some GMPEs also need as input the depth of the hypocenter or the
width of the rupture. For each earthquake scenario, these additional input parameters are
calculated for each realization of the rupture geometry and hypocenter and site locations,
and their values are then used for comparison of the synthetic motions with the GMPEs.
Moreover, the NGA-West2 GMPEs (except I14) differentiate between sites located on
the HW, that is, on the down-dip side of the top of rupture, and sites located on the foot-
wall (FW). The HW effect refers to the increase in ground motion amplitudes at short peri-
ods observed at sites located at short distances on the HW side of the rupture, as compared
to sites located at the same distance but on the FW side (Donahue and Abrahamson,
2014). This effect is accounted for in the NGA-West2 GMPEs using an input HW factor,
herein referred to as HW (= 1 for HW sites; = 0 for FW sites). The simulation procedure
of Dabaghi and Der Kiureghian does not explicitly distinguish between HW and FW sites.
Thus, the simulated motions are compared with the GMPEs for random HW–FW config-
urations. Default values are assigned for the other input parameters that remain unspeci-
fied in the various GMPEs.
Figure 11 shows the comparison for two of the considered near-fault scenarios for ran-
dom HW–FW configurations, one for an SS fault (F = 0) and one for an REV fault
(F = 1). In addition to the median and median plus/minus one standard deviation of the
synthetic ground motion spectra, we show the median spectra for the subsets of pulse-like
(P) and non-pulse-like (NP) synthetic motions. The generated synthetic ground motions
show good agreement with the NGA-West2 GMPEs. At most periods and for both the
strike-slip (F = 0) and reverse (F = 1) faulting mechanisms, the median spectrum of all the
Daoud et al. 19

Figure 12. Median of the RotD50 spectra for 600 synthetic motions versus median spectra predicted
by the NGA-West2 model (a) for scenarios having the same F, RRUP , Vs30 , and ZTOR values but different
Mw values, (b) for scenarios having the same F, Mw , Vs30 , and ZTOR values but different RRUP values, and
(c) for scenarios having the same F, Mw , RRUP , and ZTOR values but different VS30 values.

simulated motions falls within the range spanned by the median spectra of the five
GMPEs. Moreover, the median level of the simulated motions falls within the median plus
and minus one standard deviation levels predicted by the NGA-West2 model. For strike-
slip faults, the synthetic motions tend to predict larger spectral ordinates at periods
between 0.1 and 0.3 s and at longer periods (greater than 1 s) and to predict smaller spec-
tral ordinates at other periods; see Figure 11 (left). However, the differences are not large.
These observations are mostly consistent among all strike-slip scenarios. For the particular
reverse faulting scenario shown in Figure 11 (right), the synthetic motions tend to predict
larger spectral ordinates only at longer periods (greater than 1 s) and smaller spectral ordi-
nates at all other periods. While these differences are larger than the ones observed for the
strike-slip scenario, they are still not significant. Other reverse faulting scenarios exhibit
slightly different trends, whereby some have more significant differences between the simu-
lated motions and the GMPEs, while others show better agreement between the two. For
events with Mw ø 7, the simulations tend to predict spectral ordinates that are in overall
good agreement with the NGA-West2 model at longer periods, but predict smaller spec-
tral ordinates at shorter periods (less than 1 s).
It can be observed from Figure 11 that, for both fault types, the median spectrum of the
subset of simulated pulse-like motions has considerably larger amplitudes than the median
spectrum for all simulated motions as well as the median spectrum of the NGA-West2
model for periods greater than 1 s. This is due to the forward directivity effect, which is
present in the latter spectra only in a weighted average sense. On the other hand, the med-
ian spectrum of the subset of simulated non-pulse motions has smaller amplitudes than the
median spectrum of all simulated motions and smaller values than the median spectrum of
the NGA-West2 model at periods 0.5–3 s. This is due to the absence of long-period pulses
in these simulated motions. Similar trends are observed for the other scenarios.
Figure 12a to c show the effect on the median spectra of the simulated ground motions
and the NGA-West2 model of varying Mw , RRUP , and Vs30 , respectively, while keeping
other scenario parameters fixed. It is observed in Figure 12a that the effect of varying Mw
20 Earthquake Spectra 00(0)

is not captured similarly by the synthetic ground motions and the NGA-West2 model.
For the NGA-West2 model, as Mw increases, spectral amplitudes increase at all periods,
whereas for the synthetic motions spectral amplitudes increase with Mw at longer periods,
while magnitude saturation is observed at lower periods. On the contrary, Figure 12b
shows that, as RRUP decreases, the spectral amplitudes increase at all periods for both the
NGA-West2 model and the synthetic motions. Moreover, for the scenarios considered in
Figure 12b, the period-dependent differences between the NGA-West2 model and the syn-
thetic motions are similar to the differences observed in Figure 11. Finally, it is observed
in Figure 12c that the effect of varying Vs30 is not captured similarly by the NGA-West2
model and the synthetic ground motions. For the NGA-West2 model, as Vs30 decreases
spectral amplitudes increase at all periods, but more at the longer periods. For the syn-
thetic motions, the effect of Vs30 is negligible at shorter periods and spectral amplitudes
increase as Vs30 decreases only at periods longer than about 0.5 s. The trends that the syn-
thetic motions show with Vs30 at lower periods are comparable to the ones observed by
CB14. Since both the simulated motions and the NGA-West2 model are based on fitting
assumed models to data of recorded motions, it is not possible to ascertain as to which of
the two produces more accurate trends. Nevertheless, it is noted that the differences
between the two models and their trends are relatively small. And although some differ-
ences are observed in the median spectral shapes between simulated motions and the
NGA-West2 model, similarities are observed between the spectral shapes of simulated
motions and those of the ASK14 and BSSA14 GMPEs, especially at shorter periods.

Conclusion
The site-based stochastic model and simulation procedure proposed by Dabaghi and Der
Kiureghian (2014, 2017, 2018) requires information about the source, site, and source-to-
site geometry, namely the parameters F, Mw , ZTOR , RRUP , Vs30 , sor d, and uor f. The resulting
simulations account for the near-fault rupture directivity effect and capture the natural
variability of real ground motions. They also make the crucial distinction between pulse-
like and non-pulse-like ground motions. However, the directivity parameters sor d and uor f
entail information about the rupture geometry and hypocenter location that may not be
available to users of the procedure, particularly design engineers. In this article, a proce-
dure to simulate the rupture geometry and hypocenter location and thus generate synthetic
motions with random directivity effect for a given earthquake scenario is proposed.
Using data from a subset of shallow crustal earthquakes in active tectonic regions with
5:2<Mw <7:9 in the NGA West2 database, new predictive models are developed for the
rupture geometry parameters (ZTOR , LR , and WR ) as a function of F and Mw . These models
are found to be consistent with existing models. Important features of the developed mod-
els include: a two-part model consisting of a logistic regression model and a linear regres-
sion model to account for the zero-inflation in ZTOR values; account for the variability and
correlations between the rupture geometry parameters; and models for LR and WR that dis-
tinguish between buried ruptures (ZTOR .0) and surface ruptures (ZTOR = 0). These models
reveal several important observations about the scaling of rupture geometry with magni-
tude. For example, the models indicate that the scaling of the rupture width with magni-
tude is different for buried and surface ruptures, while the scaling of the rupture length
with magnitude is not affected by ZTOR .
The first simulation procedure (P1) developed in this article considers the case where
only the contributing fault, earthquake magnitude, and site location and characteristics
Daoud et al. 21

are known. This procedure is appropriate for design or assessment situations. Using
Monte Carlo simulation, random rupture geometries and hypocenter locations are gener-
ated according to their probability distributions and are used to calculate the correspond-
ing randomized rupture directivity conditions. The location of the rupture within the fault
plane is assumed to follow a uniform distribution, while the hypocenter location is
assumed to follow the distributions proposed by Mai et al. (2005). The procedure is illu-
strated by simulating near-fault ground motions at a site located in downtown Los
Angeles.
The second simulation procedure (P2) takes as input the type-of-faulting, the earth-
quake magnitude, the source-to-site distance, and the site characteristics. In this case, the
rupture geometry as well as both the hypocenter and site locations are randomized. This
procedure allows comparison with existing GMPE models. The procedure is illustrated for
several earthquake scenarios defined by the set of parameters F, Mw , ZTOR , RRUP , and Vs30 .
Results show general agreement between the two, but also some differences. For instance,
synthetic ground motions are generally consistent with the GMPEs for random HW–FW
configurations, but they are not able to capture the HW effect, though this effect is rather
small. On the contrary, the synthetic motions explicitly account for the directivity effect,
while GMPEs account for this effect in an average sense.
The proposed predictive models of rupture geometry and simulation procedure P1 are
necessary for structural engineers who want to generate synthetic motions to perform seis-
mic design or assessment studies at a near-fault site with known location and Vs30 . They
only need information about the magnitude (Mw ) and the source fault of the earthquakes
that contribute most to the hazard at their site of interest. This information can be
obtained from deaggregation of probabilistic seismic hazard analysis results.

Acknowledgments
We thank Ibrahim Alameddine and Tom Shantz for their insightful feedback and suggestions, and
Brian Chiou for providing his racetrack generation codes. We also thank the anonymous reviewers
whose constructive comments improved the paper.

Declaration of conflicting interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/
or publication of this article.

Funding
The author(s) disclosed receipt of the following financial support for the research, authorship, and/
or publication of this article: This work was supported by the State of California through the
Transportation System Research Program of the Pacific Earthquake Engineering Research Center
(PEER). Any opinions findings, and conclusion or recommendations expressed in this material are
those of the author(s) and do not necessarily reflect those of the funding agency.

ORCID iD
Mayssa Dabaghi [Link]
22 Earthquake Spectra 00(0)

Data and resources

The rupture geometry data used to develop the models in this study can be obtained from the NGA-
West2 flatfiles ([Link]

Supplemental material
Supplemental material for this article is available online.

References
Abrahamson NA, Silva WJ and Kamai R (2014) Summary of the ASK14 ground motion relation for
active crustal regions. Earthquake Spectra 30(3): 1025–1055.
Boore DM (2010) Orientation-independent, nongeometric-mean measures of seismic intensity from
two horizontal components of motion. Bulletin of the Seismological Society of America 100(4):
1830–1835.
Boore, D. M., Stewart, J. P., Seyhan, E., and Atkinson, G. M., 2014. Nga-west2 equations for
predicting pga, pgv, and 5% damped psa for shallow crustal earthquakes, Earthquake Spectra
30(3), 1057–1085.
Bozorgnia Y, Abrahamson NA, Atik LA, Ancheta TD, Atkinson GM, Baker JW, Baltay A, Boore
DM, Campbell KW, Chiou BSJ, Darragh R, Day SM, Donahue JL, Graves RW, Gregor N,
Hanks TC, Idriss IM, Kamai R, Kishida T, Kottke A, Mahin SA, Rezaeian S, Rowshandel B,
Seyhan E, Shahi SK, Shantz T, Silva W, Spudich P, Stewart JP, Watson-Lamprey J, Wooddell K
and Youngs R (2014) NGA-West2 research project. Earthquake Spectra 30(3): 973–987.
Campbell KW, Abrahamson N, Power M, Chiou BSJ, Bozorgnia Y, Shantz T and Roblee C (2009)
Next Generation Attenuation (NGA) project: Empirical ground motion prediction equations for
active tectonic regions. In: Proceedings of the 6th international conference on urban earthquake
engineering, Tokyo, Japan, 3–4 March.
Campbell KW and Bozorgnia Y (2014) NGA-West2 ground motion model for the average
horizontal components of PGA, PGV, and 5% damped linear acceleration response spectra.
Earthquake Spectra 30(3): 1087–1115.
Chiou BSJ and Youngs RR (2014) Update of the Chiou and Youngs NGA model for the average
horizontal component of peak ground motion and response spectra. Earthquake Spectra 30(3):
1117–1153.
Dabaghi M and Der Kiureghian A (2014) Stochastic modeling and simulation of near-fault ground
motions for performance-based earthquake engineering. PEER report no. 2014/20, December.
Berkeley, CA: University of California, Berkeley.
Dabaghi M and Der Kiureghian A (2017) Stochastic model for simulation of near-fault ground
motions. Earthquake Engineering & Structural Dynamics 46(6): 963–984.
Dabaghi M and Der Kiureghian A (2018) Simulation of orthogonal horizontal components of near-
fault ground motion for specified earthquake source and site characteristics. Earthquake
Engineering & Structural Dynamics 47: 1369–1393.
Dalguer LA, Miyake H, Day SM and Irikura K (2008) Surface rupturing and buried dynamic-
rupture models calibrated with statistical observations of past earthquakes. Bulletin of the
Seismological Society of America 98(3): 1147–1161.
Daoud Y, Dabaghi M and Der Kiureghian A (forthcoming) Chapter 3: Simulation of near-fault
ground motions for randomized directivity parameters. In: PEER Report: Improvements to the
Stochastic Modeling and Simulation of Ground Motions for Use in PBEE.
Donahue JL and Abrahamson NA (2014) Simulation-based hanging wall effects. Earthquake Spectra
30(3): 1269–1284.
Ellsworth WL (2003). Appendix D: Magnitude and area data for strike slip earthquakes in: US
Geological Survey Open File Report no. 03-214: Earthquake Probabilities in San Francisco Bay
Region: 2002–2031, by Working Group on California Earthquake Probabilities, Menlo Park,
CA.
Daoud et al. 23

Field EH, Dawson TE, Felzer KR, Frankel AD, Gupta V, Jordan TH, Parsons T, Petersen MD,
Stein RS, Weldon R and Wills CJ (2009) Uniform California earthquake rupture forecast, version
2 (UCERF 2). Bulletin of the Seismological Society of America 99(4): 2053–2107.
Goda K, Yasuda T, Mori N and Maruyama T (2016) New scaling relationships of earthquake
source parameters for stochastic tsunami simulation. Coastal Engineering Journal 58(3):
1650010-1–1650010-40.
Graves R and Pitarka A (2016) Kinematic ground-motion simulations on rough faults including
effects of 3D stochastic velocity perturbations. Bulletin of the Seismological Society of America
106(5): 2136–2153.
Graves R, Jordan TH, Callaghan S, Deelman E, Field E, Juve G, Kesselman C, Maechling P, Mehta
G, Milner K, Okaya D, Small P and Vahi K (2011) CyberShake: A physics-based seismic hazard
model for Southern California. Pure and Applied Geophysics 168(3–4): 367–381.
Gupta ID (2006) Defining source-to-site distances for evaluation of design earthquake ground
motion. In: Proceedings of the 13th symposium on earthquake engineering, Indian Institute of
Technology Roorkee, Roorkee, India, 18–20 December.
Hanks TC and Bakun WH (2002) A bilinear source-scaling model for M-log a observations of
continental earthquakes. Bulletin of the Seismological Society of America 92(5): 1841–1846.
Hartzell S, Guatteri M, Mai PM, Liu P-C and Fisk MR (2005) Calculation of broadband time
histories of ground motion, Part II: Kinematic and dynamic modeling using theoretical Green’s
functions and comparison with the 1994 northridge earthquake. Bulletin of the Seismological
Society of America 95(2): 614–645.
Idriss IM (2014) An NGA-West2 empirical model for estimating the horizontal spectral values
generated by shallow crustal earthquakes. Earthquake Spectra 30(3): 1155–1177.
Kagawa T, Irikura K and Somerville PG (2004) Differences in ground motion and fault rupture
process between the surface and buried rupture earthquakes. Earth, Planets and Space 56(1):
3–14.
Kaklamanos J, Baise LG and Boore DM (2011) Estimating unknown input parameters when
implementing the NGA ground-motion prediction equations in engineering practice. Earthquake
Spectra 27(4): 1219–1235.
Leonard M (2010) Earthquake fault scaling: Self-consistent relating of rupture length, width, average
displacement, and moment release. Bulletin of the Seismological Society of America 100(5A):
1971–1988.
Mai PM, Spudich P and Boatwright J (2005) Hypocenter locations in finite-source rupture models.
Bulletin of the Seismological Society of America 95(3): 965–980.
Mavroeidis GP and Papageorgiou AS (2003) A mathematical representation of near-fault ground
motions. Bulletin of the Seismological Society of America 93(3): 1099–1131.
Pacific Earthquake Engineering Research Center (PEER) (2015) NGA-West2 database. Available at:
[Link] (accessed October 2018).
Pitarka A, Graves R, Irikura K, Miyake H and Rodgers A (2017) Performance of Irikura recipe
rupture model generator in earthquake ground motion simulations with Graves and Pitarka
hybrid approach. Pure and Applied Geophysics 174(9): 3537–3555.
Rezaeian S and Der Kiureghian A (2008) A stochastic ground motion model with separable
temporal and spectral nonstationarities. Earthquake Engineering & Structural Dynamics 37(13):
1565–1584.
Rezaeian S, Petersen MD, Moschetti MP, Powers P, Harmsen SC and Frankel AD (2014)
Implementation of NGA-West2 ground motion models in the 2014 U.S. National Seismic
Hazard Maps. Earthquake Spectra 30(3): 1319–1333.
Schwartz DP (2018) Review: Past and future fault rupture lengths in seismic source
characterization—The long and short of it. Bulletin of the Seismological Society of America
108(5A): 2493–2520.
Shahi SK and Baker JW (2014) An efficient algorithm to identify strong-velocity pulses in
multicomponent ground motions. Bulletin of the Seismological Society of America 104(5):
2456–2466.
24 Earthquake Spectra 00(0)

Shapiro SS and Wilk MB (1965) An analysis of variance test for normality (complete samples)y.
Biometrika 52(3–4): 591–611.
Shaw BE (2009) Constant stress drop from small to great earthquakes in magnitude-area scaling.
Bulletin of the Seismological Society of America 99(2A): 871–875.
Somerville PG, Smith NF, Graves RW and Abrahamson NA (1997) Modification of empirical
strong ground motion attenuation relations to include the amplitude and duration effects of
rupture directivity. Seismological Research Letters 68(1): 199–222.
Southern California Earthquake Center (SCEC) (2018) CyberShake. Available at: https://
[Link]/scecpedia/CyberShake (accessed 3 March 2019).
Spudich P, Bayless JR, Baker JW, Chiou BSJ, Rowshandel B, Shahi SK and Somerville PG (2013)
Final report of the NGA-West2 directivity working group. PEER report no. 2013/09, May.
Berkeley, CA: University of California, Berkeley.
Spudich P, Rowshandel B, Shahi SK, Baker JW and Chiou BSJ (2014) Comparison of NGA-West2
directivity models. Earthquake Spectra 30(3): 1199–1221.
Thingbaijam KKS, Mai PM and Goda K (2017) New empirical earthquake source-scaling laws.
Bulletin of the Seismological Society of America 107(5): 2225–2246.
US Geological Survey (USGS) (2019) Unified hazard tool (dynamic conterminous US 2008 edition).
Available at: [Link] (accessed 3 March 2019).
Watts AB and Burov EB (2003) Lithospheric strength and its relationship to the elastic and
seismogenic layer thickness. Earth and Planetary Science Letters 213(1–2): 113–131.
Wells DL and Coppersmith KJ (1994) New empirical relationships among magnitude, rupture
length, rupture width, rupture area, and surface displacement. Bulletin of the Seismological
Society of America 84(4): 974–1002.
Witte RS and Witte JS (2016) Statistics. Hoboken, NJ: Wiley.