J Seismol
DOI 10.1007/s10950-008-9090-6

ORIGINAL ARTICLE

Broadband strong motion simulation in layered half-space
using stochastic Green’s function technique
Y. Hisada

Received: 30 March 2007 / Accepted: 7 January 2008
# Springer Science + Business Media B.V. 2008

Abstract The stochastic Green’s function method,which simulates one component of the far-field S-waves
from an extended fault plane at high frequencies(Kamae et al., J Struct Constr Eng Trans AIJ, 430:
1–9, 1991), is extended to simulate the three components of the full waveform in layered half-spaces forbroadband frequency range. The method firstly computes ground motions from small earthquakes, whichcorrespond to the ruptures of sub-faults on a fault planeof a large earthquake, and secondly constructs thestrong motions of the large earthquake by superposing
the small ground motions using the empirical Green’s
function technique (e.g., Irikura, Proc 7th Japan
Earthq Eng Symp, 151–156, 1986). The broadband
stochastic omega-square model is proposed as themoment rate functions of the small earthquakes, inwhich random and zero phases are used at higher andlower frequencies, respectively. The zero phases are
introduced to simulate a smooth ramp function of the
moment function with the duration of 1/fc s (fc: thecorner frequency) and to reproduce coherent strong
motions at low frequencies (i.e., the directivity pulse).
As for the radiation coefficients, the theoretical valuesof double couple sources for lower frequencies and thetheoretical isotropic values for the P-, SV-, and SH-

Y. Hisada (*)
Department of Architecture, Kogakuin University,
Nishi-Shinjuku 1-24-2, Shinjuku,
Tokyo 163-8677, Japane-mail: hisada@cc.kogakuin.ac.jp

waves (Onishi and Horike, J Struct Constr Eng Trans
AIJ, 586:37–44, 2004) for high frequencies are used.
The proposed method uses the theoretical Green’s
functions of layered half-spaces instead of the far-fieldS-waves, which reproduce the complete waves including the direct and reflected P- and S-waves and surface
waves at broadband frequencies. Finally, the proposed
method is applied to the 1994 Northridge earthquake,and results show excellent agreement with the observation records at broadband frequencies. At the sametime, the method still needs improvements especiallybecause it underestimates the high-frequency verticalcomponents in the near fault range. Nonetheless, themethod will be useful for modeling high frequency
contributions in the hybrid methods, which use
stochastic and deterministic methods for high and low
frequencies, respectively (e.g., the stochastic Green’s
function method+finite difference methods; Kamae
et al., Bull Seism Soc Am, 88:357–367, 1998; Pitarka
et al., Bull Seism Soc Am 90:566–586, 2000), because
it reproduces the full waveforms in layered media
including not only random characteristics at higher
frequencies but also theoretical and deterministic
coherencies at lower frequencies.
Keywords Broadband strong motion simulation .
Omega-squared model .
Green’s function of layered half-spaces .
Stochastic Green’s function method .
Empirical Green’s function method .
The scaling law . 1994 Northridge earthquake

DO09090; No of Pages

J Seismol

1 Introduction
Hybrid methods that combine deterministic andstochastic estimations for lower and higher frequencies are widely used to simulate broadband strong
ground motions (e.g., Kamae et al. 1998; Pitarka et al.
2000). The matching frequencies between the two
estimations are usually around 0.5–2 Hz, which are
the most important frequency range for engineering
structures. However, the two estimations sometimes
show large differences in amplitudes, arrival times,
and waveforms in that frequency range because ofcompletely different methodologies. Two approaches
have been proposed to remedy this problem andconnect the two results more smoothly. One approachis to introduce stochastic randomness in deterministic
methods at higher frequencies. For example, Herrero
and Bernard (1994) proposed the k-squared model,
which is an omega-squared model based on a kinematic fault model using randomness in slip distribution
at higher frequencies. Hisada (2000, 2001) modified
the model by introducing pseudo-dynamic slip andthe k-squared random distributions not only in slip but
also in rupture time at higher frequencies. Morerecently, various similar models considering pseudodynamic slips and randomness for slip and rupture
time (or rupture velocity) have been proposed (e.g.,
Guatteri et al. 2004; Hartzell et al. 2005). These
approaches are theoretically sound, but rather time
consuming and not efficient at very high frequencies,
especially when considering realistic layered halfspaces, because huge numbers of Green’s functions
are necessary to maintain the continuity of the rupture
front (about 8–10 source points per wavelength).
The second approach introduces deterministic coherency in stochastic methods at lower frequencies.
Kamae et al. (1991) proposed the stochastic Green’s
function method, in which the far-field S-waves fromthe stochastic point source (Boore 1983) are superposed on an extended fault plane using the empiricalGreen’s function technique (Irikura 1986). To extend
this method to lower frequencies, Kagawa (2004)
proposed a trial-and-error technique, in which several
stochastic Green’s functions are simulated, and one canchoose synthetics that are consistent with deterministic
results. Onishi and Horike (2004) improved the methodby introducing the theoretical radiation coefficients of
the P-, SV-, and SH-waves and a ray tracing techniquein layered half-spaces. These stochastic approaches aremuch faster to compute than the theoretical approaches
because it basically requires only one source point per
sub-fault using simple Green’s functions.
The purpose of this paper is to improve the stochastic
Green’s function method by introducing deterministic
coherencies at lower frequencies and by computing thecomplete Green’s functions of layered half-spaces. First,
we propose the broadband stochastic omega-squaremodel, whose phases are random at higher frequencies
and zero at lower frequencies, to reproduce both therandom characteristics at higher frequencies and thedeterministic coherencies (e.g., the directivity pulses) atlower frequencies. Second, we introduce the theoreticalGreen’s function of layered half-spaces using the wavenumber integration technique based on the reflection/transmission (R/T) matrix (Hisada 2001), which isnumerically stable up to very high frequencies. Third,
we propose a method to compute the theoretical P-,
SV-, and SH-radiations at high frequencies (Onishi andHorike 2004), by introducing P-, SV-, and SH-sources.
Finally, we check the validity of the method bycomparing observations with results from simulations
for the 1994 Northridge earthquake.

2 Formulation of the broadband stochastic Green’s
function method in layered half-space
2.1 Formulation of the stochastic Green’s function
methodWe briefly explain the original stochastic Green’s
function method (Kamae et al. 1991; see also Kamae
et al. 1998). In the method, one first generates groundmotions from small earthquakes, which correspond to
sub-faults of a large earthquake as shown in Fig. 1.
Next, one constructs the strong motions of the large
earthquake by superimposing those of the small
earthquakes, using the empirical Green’s function
technique (Irikura 1986; Yokoi and Irikura 1991).
Although the original formulation is written in thetime domain, we will express it in the frequency
domain, as follows:
UkL ð5 Þ ¼

N X
NXξ¼1 η¼1

SF ð5 Þ  Cξη  Ukξη
ð5 Þ  exp i5  tξη ð1Þ

ðk ¼ x; y; zÞ;J Seismol

as semi-empirical Green’s functions, whereas Kamae
et al. (1991) used the stochastic Green’s function of
the far field S-waves in a homogeneous full-space(Boore 1983). In addition, those methods assume
random phases for the source spectra; this restricts theextension to lower frequencies. We extend the methodby introducing coherencies in source spectra at lower
frequencies and by using the theoretical Green’s
function in a layered half-space.
Using the representation theorem of a point
dislocation source (e.g., Aki and Richards 1980), thevelocities at an observation point Y from the small
earthquakes are represented as follows:::
UkS ðXO ; 5 Þ ¼ M S ð5 Þ ei nj þ ej ni Uik;j ðXO ; XS ; 5 Þ ð4Þ

Fig. 1 Fault parameters of a large fault and sub-faults

where ω is the circular frequency, and the superscripts
L and S stand for the large and small earthquakes,
respectively. The subscript k corresponds to a component of the x-, y-, or z-axes. ξ and η are the sub-fault
numbers along the length and the width, and txm is the
rupture time of the ξη sub-fault (see Fig. 1). The
function F is introduced to modify the slip function of
the small earthquakes into that of the large earthquake, as we will explain in Section 2.4. N is thenumber of sub-faults along the length and the widthand is also the scaling parameter between the large
and small earthquakes. This method assumes thefollowing scaling relation between the small and large
earthquakes (Irikura 1983; Yokoi and Irikura 1991).
LL W L DL τ L¼ S¼ S¼ S¼
LS
W
D
τ



M0L
CM0S

1=3
N

ð2Þ

where, L, W, D, and τ are the length, the width, theslip, and the rise time of the faults, respectively. Cxm
and C in Eqs. 1 and 2 are the ratio of the stress drops
between the large and small earthquakes (Yokoi and
Irikura 1991).
C¼

$σL
$σS

ð3Þ

2.2 Broadband omega-square model for small
earthquakes and moment rate function
Regarding the Green’s functions from small earthquakes (UkS ), Irikura (1986) used observation records

ðk ¼ x; y; zÞwhere, XO and XS are observation and source points,
respectively (see Fig. 1). The subscripts i, j, and kcorrespond to components of the x, y, or z-axes, where
the summation convention is used. ei and ni are the ith
component of the unit vectors in the fault slip and thefault normal direction, respectively. Uik,j is the derivative of the Green’s function with respect to the jth
direction. The theoretical Green’s functions of thelayered half-space are formulated by the wavenumber integration method based on the R/T matrix
method, and easily computed even at extremely highfrequencies(e.g., Hisada 1995).
:
M S in Eq. 4 is the source spectrum of a small
earthquake, which is the Fourier spectrum of the
moment-rate function. Assuming the omega-squaredmodel (Brune 1970; Boore 1983), its amplitude
spectrum at a frequency, f, is expressed as follows: :S 
 M ð f Þ ¼

 S 
M0S
  S 2 P f ; fmax
1 þ f fC

ð5Þ

where, M0S is the seismic moment of the small
earthquake (in dyne-cm), and is expressed as follows
using Eq. 2.

M0S  M0L CN 3 ¼ μDL LL W L CN 3

ð6Þ

fCS in Eq. 5 is the corner frequency (Hz), expressed asfollows (Brune 1970), 1=3fCS ¼ 4:9  106 VS $σS M0S

ð7Þ

J Seismol

1
P ðf ; fmax Þ ¼ pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ
1 þ ðf =fmax Þ 2n

ð8Þ

3.5

fc=1 Hz

Duration= 1/ fc s

moment-rate function

where, VS is the shear wave velocity of the source
layer (in km/s), and ΔσS is the stress drop (inbar=0.1 MPa). The function P in Eq. 5 is a high-cut
Sfilter due to fmax,

Moment-Rate Function (original)

3

Moment-Rate Function (0.5/fc s delay)
2.5
2Moment Function (original)
Moment Function (0.5/fc s delay)

1.5
1
0.5
0

where, n=4 is recommended by Boore (1983). On theother hand, Tsuruki et al. (2006) proposed fmax =6.0
Hz and n=1.55 for crustal earthquakes, using various
observation data in Japan.
We propose the broadband omega-square modelusing not only random phases from −π to +π at higher
frequencies but also zero phases at lower frequencies.
The use of zero phases is one of the simplest
techniques to reproduce a coherent moment-ratefunction and thus to simulate the coherent waves at
low frequencies. To explain this, we first derive thetheoretical moment function of the omega-squaremodel by assuming zero phases. Equation 5 with zero
phases is transformed into the following moment-ratefunction, using the Fourier inverse-transform.
:
M S ðt Þ ¼ π  M0S fCS exp 2π  fCS jt j

ð9Þ

where we assume that P=1 for simplicity; numericalresults show that the inclusion of P gives a smoothershape of the moment-rate function. The moment
function corresponding to Eq. 9 is expressed asfollows:(
M ðt Þ ¼
S
0:5M0S exp þ2π  fCS t ; ðt  0Þ
0:5M0S 2  exp 2π  fCS t ; ðt  0Þ
ð10Þ

Figure 2 shows an example of the moment-ratefunction and the corresponding moment functions forthe case of fCS ¼ 1 Hz and M0S ¼ 1. The original
moment-rate function (the thin blue line with “original”
in the figure) is a simple axial symmetry function
peaked at t=0 (s), and the moment function (thick blueline) is a smoothed ramp function with duration
 S(risetime) of 1 fCS (s). Since they start at about
0:5
fC (s),
 S
we introduce a delay time of þ0:5 fC (s) for actual

-1

-0.5

0

0.5
time (s)

1

1.5

2

Fig. 2 Normalized moment-rate and moment functions using
the omega-square model with zero phases

applications, as shown by the red lines with “0.5/fc sdelay” in the figure.
Next, we show an example of the broadband
moment-rate and moment functions using both zeroand random phases at lower and higher frequencies,
respectively. Figure 3 shows results of the case forfCS ¼ 0:25 Hz and M0S ¼ 1. Figure 3a is a broadband
phase spectrum, in which we use random phases
between −π and +π at frequencies higher than fr (=2Hz), and the phases gradually tend to zero at
frequencies lower than fr. Figure 3b shows the Fourier
amplitudes of the source spectra, in which the blackline is the theoretical omega-squared model of Eq. 5,
and the red line represents the iterated results obtainedby using the phase spectrum of Fig. 3a. Figure 3cshows the moment-rate functions with a duration of
about 1 fCS ¼ 4 s, which correspond to the displacement waveforms in far-field. The black and red lines
represent the results using the zero phases and thephases of Fig. 3a, respectively. In the latter result, weclip the amplitudes before t=0 (s) and negativevalues. Figure 3d shows the moment-functions, andFig. 3e and f are the first- and second-orderderivatives of the moment-rate functions, whichcorrespond to the velocities and the accelerations in
the far-field, respectively. The results without random
phases show unrealistic waveforms, especially in theacceleration. On the other hand, the broadband
omega-square model reproduces not only the coherent
moment-rate and moment functions at lower frequencies (see Fig. 3c,d) but also the random waveforms at
high frequencies (see Fig. 3e,f).

J Seismol
Fig. 3 Moment andmoment-rate functions
using the broadband stochastic source model(fc=0.25 Hz). a Fourierphase spectrum (fr =2 Hz).
b Omega-squared Fourier
amplitude spectra.
c Moment-rate function(far-field displacement).
d Normalized moment
function. e First derivativeof moment-rate function(far-field velocity). f Second
derivative of moment-ratefunction (far-field
acceleration)

2.3 High-frequency source model for small
earthquakes
In this section, we will describe the procedure to
introduce the radiation coefficients of P-, SV-, andSH-waves in a layered half-space at high frequencies
and to combine these waves with the theoreticalradiation pattern of a double couple source at low
frequencies. In the case of the far-field in homoge-

neous media, it is easy to combine the radiation
coefficients at high frequencies with the theoreticalradiation pattern at low frequencies because the P- andS-waves are expressed independently (e.g., Pitarka et al.
2000).
For the high frequency radiation coefficients,
Boore (1983) used the homogenousradiation coeffipﬃﬃﬃﬃﬃﬃﬃﬃ
cient of the S-wave (¼ 2=5, Prob. 4.6 in Aki andRichards 1980), by averaging the far-field S-wave

J Seismol

over whole focal sphere in a homogenous space.
Boore and Boatwright (1984) calculated numericallycoefficients of P-, SV-, and SH-waves for various
combinations of source and takeoff angles. Recently,
Onishi and Horike (2004) derived the analyticalradiation coefficients of the P-, SV-, and SH-waves
in a homogenous space as follows:

1 2
cos2 1 
2 3

2
þ cos2 δ
51
14
sin2 1 þ sin2 ð2δÞ
3
5

1 2
cos2 1  1 þ sin2 δ
FSH ¼
2 3
þ

þ

:
UkSS ðXO ; 5 Þ ¼ :
FS  M S ð5 Þ ei nj þ ej ni Uik;j ðXO ; XS ; 5 Þð13Þ1
FP ¼ ð4=15Þ =2
FSV ¼

be perpendicular to the lines from the source tothe observation points. For these cases, the highfrequency waves from the S sources are expressedby the following equations.

1
sin2 1  1 þ cos2 ð2δÞ
3

1=
2

ð11Þ

1=
2

where, δ and 1 are the dip and rake angles of a source.
FP is the same as Prob. 4.6 in Aki and Richards(1980). It is easily confirmed that the root meansquare value of FSV andpFﬃﬃﬃﬃﬃﬃﬃ
SH ﬃ is equal to the S-wave
radiation coefficient (¼ 2=5).
To apply these coefficients to a layered halfspace, we need to generate independently the P-,
SV-, and SH-waves. We compute them approximately by modifying the representation theorem of
Eq. 4. First, we employ the explosive source forgenerating the P-wave, using the following equation(see Fig. 4a):::
UkSP ðXO ; 5 Þ ¼ FP  M S ð5 Þ  Uik;i ðXO ; XS ; 5 Þð12ÞNext, we use dip and strike slip sources togenerate the SH- and SV-waves, respectively. As
shown in Fig. 4b and c, their fault planes are set toFig. 4 The explosive, dip,and strike slip sources togenerate the high-frequency
P-, SV- and SH-waves,respectively. a Explosive
source. b Dip slip source.
c Strike slip source

where, FS is the radiation coefficient of the SV- andSH-waves in Eq. 11, and the slip vectors (ei) are thoseof the dip and strike slip sources shown in the Fig. 4b
and c for the SV- and SH-waves, respectively.
It is well know that two horizontal components do
not show large differences in amplitudes at highfrequencies as compared with those at low frequencies. To make the two horizontal components nearly
equal at high frequencies, first, we divide the radial
components of the results from the P- and SV-sourcespﬃﬃﬃ
equally into the X- and Y-components by
2.
Similarly, we divide the transverse components fromthe SH-source into the two components equally. We
generate two sets of the broadband omega-squaremodel with different random phases and use theirphases in the two components such that they are
independent of each other. Finally, we superpose
those high-frequency results with the low frequency
results from Eq. 4, by using high-pass and low-passfilters, respectively; we use a trapezoidal filter withtapering amplitudes from one to zero at matchingfrequencies.
2.4 F-functions and scaling laws of source spectra
It only remains to describe the procedure for selecting
the F-function in Eq. 1 to ensure that the scaling lawof the source spectra between the small and large

J Seismol

earthquakes will be obeyed. Under the assumptions ofpoint sources and the omega-square model withsimilar values of fmax, the scaling law is expressed
as follows (e.g., Yokoi and Irikura 1991).

8 L
M
:L
>
<CM0 S  N 3
M ðfÞ
0
:S

 2

2=

1=
M
>
C M ð f Þ : 0L fCL ¼ M0L $σL M0S 3 ¼ M0S 3  N
CM0S

fCS

CM0S

M0L $σS

CM0L

The above scaling law is easily confirmed using Eqs.
2 and 5 for low frequencies and Eqs. 2, 3, 5, and 7 forhigh frequencies.
We test two types of the F-functions.

sin ð5 τ L =2βÞi5 τ LF ð 5 Þ ¼ 1 þ ð N  1Þexp
5 τ L =2β2βð15Þ1  i5 τ S α
F ð5 Þ ¼ N 
1  i5 τ L =α
¼ 1 þ ð N  1Þ1
1  i5 τ L =αð16ÞEquation 15 is used in the semi-empirical Green’s
function method by Irikura (1986), which is the
combination of the delta function and the box-car
function with a duration of τL. On the other hand, Eq.
16 is the direct ratio of the slip functions between thesmall and large earthquakes (Onishi and Horike 2004),
in which the exponential function is assumed[f ðt Þ ¼ expðt=τ Þ=τ in the time domain, and F ðωÞ ¼
1=ð1  iωτ Þ in the frequency domain; Ben-Menahem
and Toksoz 1963]. Note that F ≈ N for low frequencies(f → 0) and F≈1 for high frequencies (for f → ∞).
We confirm that Eq. 1, with the use of the two trial
F-functions, follows the scaling law of Eq. 14 byusing Eqs. 1 and 4.
:
UL
M Lð f Þ:
 kS
C M S ð f Þ CUk
N P
N


P
S
F ð5 Þ  Cξη  Ukξη
 exp i5  tξηξ¼1 η¼1
¼
CUkS


N X
N
X



F ð5 Þ  exp i5  tξηξ¼1 η¼1

ð17Þ

ðf ! 0Þðf ! 1Þ

ð14Þ

where we assume that the focal mechanisms, Green’s
functions, and the ratios of the stress drops among thesmall earthquakes are similar. For low frequencies,
Eq. 17 is of order N3 because of the coherent
summations of N with respect to ξ (length), η (width),and F (slip ≈ N). For high frequencies, it becomes oforder N because of the random summation of N withrespect to ξ and η with F ≈ 1.
We introduce the free parameters, α and β, in
Eqs. 15 and 16 to adjust the durations of the slip
functions of the large earthquake. Figure 5a and bshows examples of the Fourier amplitudes of Ffunctions for the cases of N=10, τL =1 s, and τL =10 s.
The thin black lines and thick red lines correspond to
Eqs. 15 and 16, respectively. All the functions
become F=10 (=N) and F=1 at lower and higher
frequencies, respectively. Equation 15 shows artificial oscillations due to the presence of the sinefunction, whereas Eq. 16 shows smooth functions,which indicates that Eq. 16 is more appropriate forpractical use. The figures show F-functions for thecases of constant α and β (2β=α=1.0) and theconstant τ L =2β and τ L =α (τ L =2β ¼ τ L =α ¼ 0:5).
Frequencies between 0.5 and 2 Hz are generally
considered to correspond to the boundary between
coherent and random phases, and thus, it is probably
true that the F-functions reduce their amplitudes
from N to 1 within these frequencies. When we use
2β=α=1.0, the F-functions reduce their amplitudes
at very low frequencies especially for a large τL (=10 s,
see Fig. 5b).
Figure 6 shows various F-functions of Eq. 16,
corresponding to different values of τ L =α. Values
from 0.2 to 1.0 of this parameter seem appropriate
because they reduce the amplitudes between 0.5 and 2
Hz. Therefore, we recommend the following values
for Eqs. 15 and 16.
τ L 2β ¼ τ L α  0:2  1:0

ð18Þ

J Seismol
Fig. 5 Examples of
F-functions using Eqs. 15
and 16 (N=10). a Case forτL =1 s. b Case for τL =10 s

3 Application to strong ground motions of the 1994
Northridge earthquake
3.1 Ground motions of a small earthquake
We apply the proposed method to observed records
during the 1994 Northridge earthquake. Before
simulating the main motions, we will check groundmotions from a small earthquake using various
Green’s functions. Table 1 shows the material
properties for rock (R) and sediment (S) site models(Wald et al. 1996). Table 2 shows the source
parameters of a trial small earthquake, and Fig. 7
shows the locations of the source and the threestations, whose epicentral distances are 1, 10, and100 km along the y-axis. Figure 8 shows the fourmodels whose layered structures correspond to therock layers of Table 1. Model 1 computes the far-field12

τ = 10 s
0.5 HzTable 1 Material properties for the rock and sediment models(Wald et al. 1996)

2 Hz

F-Function

10

Layer Density(kg/m3)8

τ L /α =0.1
=0.2
=0.5
=1.0
=2.0
=5.0

6

1
2
3
4
5
6
7
8

4
2

frequency (Hz)
0
0.001

0.01

0.1

Fig. 6 F-functions using Eq. 16

S-waves at a depth of 4 km under the stations andmultiplies the 1-D amplification factor of the S-waves
using the overlying three layers. Model 2 computes
the P and S-waves in the homogenous full space
including the intermediate and near-field terms andmultiplies by the 1-D amplification factors of the Pand S-waves for the vertical and horizontal components, respectively. Models 3 and 4 compute thecomplete waves in the layered half-spaces; Model 3
uses the top four rock layers of Table 1, whereas
Model 4 uses all the layers to reproduce the reflection
waves from the Conrad and Moho discontinuities.
Figure 9a and b shows the three components of thecomputed displacements and accelerations, respectively. We use fmax =6 Hz and n=1.55 in Eq. 8, theboundary frequency (fr) between the coherent andrandom phases is 2 Hz (see Fig. 3a), and the matchingfrequencies between the theoretical and homogeneousradiation coefficients are 0.5–2.0 Hz. All the waves

1

10

1,700
1,800
2,100
2,400
2,700
2,800
2,900
3,300

Vp
Qp
(m/s)

800
1,200
1,900
4,000
5,500
6,300
6,800
7,800

40
60
100
200
400
400
600
600

Qs
Vs
(m/s)

300
500
1,000
2,000
3,200
3,600
3,900
4,500

20
30
50
100
200
200
300
300

Thickness Depth
(km)
(km)
R

S

–
–
0.5
1.0
2.5
23.0
13.0
–

0.1
0.2
0.2
1.0
2.5
23.0
13.0
–

0.1
0.3
0.5
1.5
4.0
17.0
40.0

All the Q values are assumed Qð f Þ ¼ Q0  f 0:8 , and Q( f )>Q0,
where Q0 corresponds to QP and QS in the above table

J Seismol
Table 2 Source parameters of a small earthquake
Depth

Strike Dip Rake Δσ

10 km 122°

40°

100°

fC

M0

MW

100 bar 2 Hz 6.9·1022
4.5
dyne-cmare computed up to 12.5 Hz and high-cut filtered from10 to 12.5 Hz. Note that the baselines of the waves for
Models 1–3 are shifted upward in the figures.
At station 1, which is nearly above the source, thefour models show similar results in the horizontal
components of motion because the 1-D S-waves are
dominant. However, Models 1 and 2 overestimate thehorizontal components compared with Models 3 and4 because they neglect the geometrical attenuations
within the upper three layers with the 4-km thickness.
In addition, Model 1 shows much smaller amplitudes
in the vertical components and smaller horizontal
components at lower frequencies (see the EW
components of the displacements around 4 s) thanthose of other models. This is because Model 1 lacks
the P-waves and the near-field and intermediate termsof the S-waves. On the other hand, the verticalaccelerations of all the models at Station 1 are muchsmaller than those of the horizontals because the SV
and SH sources are dominant in horizontal components than the vertical ones for stations above thesource (see Fig. 4b,c). These results appear tocontradict those of the empirical relations, in which
the vertical accelerations are generally around one
half to one third those of the horizontals. Thus, togenerate realistic vertical high frequency S-waves at a
station above a source, we may need to introduce
some modifications in the source models by considering the empirical relations in the future.

Fig. 7 Locations of thesmall earthquake (source
point) and the observation
stations

As the stations go farther away from the source(Stations 2 and 3), Models 1 and 2 show unrealistic
waveforms with much larger amplitudes and shorter
durations than those of Models 3 and 4. In particular,
they tend to overestimate the vertical components
because the SV waves in a simple homogeneous
media become dominant in vertical components at alarge distance. On the other hand, the results between
Models 3 and 4 are almost identical at Stations 1 and
2, whereas Model 4 generates more complicatedaccelerations at Station 4 than those of Model 3. Thisis because the reflected waves are generated from theConrad and Moho discontinuities (e.g., Burger et al.
1987; see accelerations for t>30 s).
3.2 Strong ground motion of the 1994 Northridge
earthquake
Finally, we simulate the strong motions during the 1994
Northridge earthquake. Figure 10a shows the locationsof the epicenter (the star in the figure), the fault plane,
which is dipping toward the southwest, and theobservation stations (triangle marks). The white andshaded areas of the figure indicate the rock and sediment areas, respectively, whose structures are shown in
Table 1 (Wald et al. 1996). Figure 10b shows the final
slip distribution of the strong-motion inversion modelby Wald et al. (1996).
We test three models; the first model is identicalwith the strong-motion inversion model by Wald et al.
(1996), where the dislocation time history for eachsub-fault is represented by the integral of threeisosceles triangles with a duration of 0.6 s and aninterval time of 0.4 s. Each triangle is allowed to havedifferent rake angles. We simulate these strong
motions up to 1.5 Hz using the layered structures
shown in Table 1. The other two models are the sameJ Seismol
Fig. 8 Four models (thenumbers in the layers correspond to those in Table 1).
Model 1 Far-field S waves +1D amplification). Model 2
P&S waves in full-space +1D amplification). Model 3
Complete waves in layered
half-space. Model 4 Complete waves in layered
half-spaceas Models 2 and 4 in Fig. 8, where we simulate strong
motions up to 12.5 Hz and use the high-cut filter from10 to 12.5 Hz. The followings are the common source
parameters in the last two models:– Scaling parameter: N=14, rupture velocity: Vr =3
km/s– Ratio of the stress drop between the small andlarge earthquakes: Chx ¼ 1

Fig. 9 Simulated waves
from the small earthquake
using the four models of
Fig. 8. Note that the baselines of the waves for
Models 1–3 are shifted
upward. a Displacement
waves. b Acceleration
waves

–
–
–
–

Seismic moments of the. small earthquakes inS¼ μDLiξη LW N 3
Eq. 6: M0ξηS
Corner frequency of the small earthquake: fCxh
inS
Eq. 7 using M0xh
High-cut filter due to fmax: fmax =6 Hz, and
n=1.55 in Eq. 8
Boundary frequency between the coherent andrandom phases: fr Hz (see Fig. 3a)

J Seismol
Fig. 9 (continued)––

Matching frequencies from the theoretical toisotropic radiation coefficients: 0.5–2.0 Hz
F-function and the adjusted rise time: Eq. 16 andτ L =α ¼ 0:5.

We use the distribution of the final slips and therake angles of Wald et al. (1996), shown in Fig. 10b.
Fig. 10 1994 Northridge
earthquake model andstrong motion stations(Wald et al. 1996). a Locations of fault, epicenter, andstations. b Slip distribution
and asperities A1 and A2

For the stress drop, we use the average values overthe fault plane, simply assuming a circular crackmodel (Eshelby 1957).
$σL ¼

7π3=2 M0L
16S 3=2

ð19Þ

J Seismol
Fig. 11 Observed and simulated waves for the 1994
Northridge earthquake. FN
and FP represent the fault
normal and parallel components (N32E and N58W),
respectively. a Observed
and simulated velocities at
U56, SYL, U03, and U53.
b Observed and simulated
accelerations at U56, SYL,
U03, and U53

where, S is the fault area. We obtain Δσ=2.98 MPa(29.8 bar) using S=L×W m2 (L=18,000 m, W=
24,000 m), and M0 =1.10×1019 N·m.
Due to limited space, we will show results at four
stations, U56, SYL, U03, and U53 (see Fig. 10a).
Figure 11a and b shows the fault normal (FN, N32E),
fault parallel (FP, N58W), and UD components of thevelocities and accelerations, respectively, at the four
stations. In each case, we plot the observed records(shown as “Observation” in the figures), the simulated
results using the original model (“Original”; Waldet al. 1996), and Models 2 and 4 (see Fig. 8). Notethat the baselines of the waves for Observation,
Original, and Model 2 are shifted upward in the
figures.
Because the original model lacks frequencies
higher than 1.5 Hz, it reproduces the observed
velocities well but not the accelerations. On the other
hand, Models 2 and 4 simulate well the horizontal
velocities and accelerations of the observations. As
for the vertical components, Model 2 generally
overestimates the velocities because the SV-waves ofModel 2 become dominant in the vertical direction ata certain distance, as seen in Section 3.1 previously.

J Seismol
Fig. 11 (continued)

On the contrary, Model 4 tends to underestimate thevertical accelerations because the SV- and SH-waves
are dominant on the horizontal components at thestations above sources.
Because U56 and SYL are located at the north ofthe epicenter, the FN components of their velocities
show the forward directivity pulses, namely, large
amplitudes and short durations. These coherent waves
are successfully reproduced by all the models. On theother hand, the velocity records of U03 and U53 near
the epicenter show near random characteristics,
namely, smaller amplitudes and longer durations.

These waves are also successfully simulated by all
the models, although the simulations generally show
shorter durations than the observations because they
do not include the scattering effects.
Figure 12 shows the Fourier amplitudes of theacceleration records (thick gray lines), the originalmodel (thick black lines), and Models 2 and 4 (thinred and blue lines, respectively). The original model
agrees well with the observations at frequencies lessthan 1 Hz, whereas Models 2 and 4 successfully
reproduce the observations at broadband frequencies.
In particular, the results of Model 4 are almostJ Seismol
Fig. 12 Observed and simulated Fourier amplitudes of
the accelerations at U56,
SYL, U03, and U53 for the
1994 Northridge earthquake

identical with those of the original model at frequencies lower than 0.2 Hz, whereas Model 2 is not able toreproduce them (e.g., the FP components at U56 and
SYL).

4 Conclusions
We propose a practical method for simulating thethree components of broadband P, S, and surface
waves in layered half-spaces from an extended fault
model, which is an extension of the original stochastic
Green’s function method (e.g., Kamae et al. 1991).
The simulated examples show that the far-field S-

waves, which are used in the original method, haveimitations not only in the near-field but also in the farfield because they exclude the near and intermediate
terms, the surface waves, and the reflected waves
from the Moho and Conrad discontinues. We also
confirm that the method successfully simulates theobserved strong motions of the 1994 Northridge
earthquake. On the other hand, we find that themethod needs improvements; first, it does not
reproduce large amplitudes in vertical high frequency
S-waves in the epicentral area because the upwardpropagating SV waves are dominant on the horizontal
component. The assumption of horizontally flatlayered structures may be an oversimplification inJ Seismol

some cases, and we may need to introduce anempirical relation between horizontal and verticalamplitudes. Second, the durations of the simulated
waves are generally shorter than those of theobservations, although we introduced the complete
waves in layered structures. To remedy this problem,
we plan to test models in the future that include thinlayers to reproduce scattering effects and/or introduce
an empirical duration function. Despite these shortcomings, the proposed method will be a more
effective tool than the original method for modelingthe high frequency contributions in the hybrid
methods (e.g., Kamae et al. 1998, Pitarka et al.
2000) because the method reproduces not only
stochastic randomness at high frequencies but also
deterministic coherencies at low frequencies, using
the complete waves in layered half-spaces.
Acknowledgements This research is partly supported by
Earthquake and Environmental Research Center (EEC) of
Kogakuin University funded by the Frontier Research Promotion Program of the Ministry of Education, Culture, Sports,
Science, and Technology of Japan (MEXT) and by the Japan
Nuclear Energy Safety Organization (JNES). The manuscript
was greatly improved by reviews from Jacobo Bielak, Kim
Olsen, and two anonymous reviewers.

References
Aki K, Richards PG (1980) Quantitative seismology. Freeman,
New York
Ben-Menahem A, Toksoz MN (1963) Source mechanism fromspectrums of long period surface waves. J Geophys Res
68:5207–5222
Boore DM (1983) Stochastic simulation of high-frequency
ground motions based on seismological models of theradiated spectra. Bull Seism Soc Am 73:1865–1894
Boore DM, Boatwright J (1984) Average body-wave radiation
coefficient. Bull Seism Soc Am 74:1615–1621
Brune JN (1970) Tectonic stress and the spectra of seismic shear
waves from earthquake. J Geophys Res 75:4997–5009
Burger RW, Somerville PG et al (1987) The effect of crustal
structure on strong ground motion attenuation relations ineastern north America. Bull Seism Soc Am 77:420–439
Eshelby JD (1957) The determination of the elastic field of anellipsoidal inclusion, and related problems. Proc Roy Soc
A241:376–396
Guatteri M, Mai P, Beroza G (2004) A pseudo-dynamic
approximation to dynamic rupture models for strong

ground motion prediction. Bull Seism Soc Am 94:2051–
2063
Hartzell S, Mariagiovanna G et al (2005) Calculation ofbroadband time histories of ground motion, Part II:
kinematic and dynamic modeling using theoretical Green’s
functions and comparison with the 1994 Northridge
earthquake. Bull Seism Soc Am 95:614–645
Herrero A, Bernard P (1994) A kinematic self-similar rupture
process for earthquake. Bull Seism Soc Am 84:1216–1228
Hisada Y (1995) An efficient method for computing Green’s
functions for a layered half-space with sources andreceivers at close repths (Part 2). Bull Seism Soc Am
85:1080–1093
Hisada Y (2000) A theoretical omega-squared model considering the spatial variation in slip and rupture velocity. Bull
Seism Soc Am 90:387–400
Hisada Y (2001) A theoretical omega-squared model considering the spatial variation in slip and rupture velocity. Part 2:
case for two-dimensional source model. Bull Seism Soc
Am 91:651–666
Irikura K (1983) Semi-empirical estimation of strong groundmotions during large earthquakes. Bull Disas Prev Res Inst
33:63–104
Irikura K (1986) Prediction of strong ground acceleration
motions using empirical Green’s function. Proc 7th Japan
Earthq Eng Symp, pp 151–156
Kagawa T (2004) Developing a stochastic Green’s function
method having more accuracy in long period range to be
used in the hybrid method. J Jpn Assoc Earthqu Eng 4:
21–32 (in Japanese)
Kamae K, Irikura K, Fukuchi Y (1991) Prediction of strong
ground motion based on scaling law of earthquake bystochastic synthesis method. J Struct Constr Eng Trans AIJ
No430:1–9 (in Japanese)
Kamae K, Irikura K, Pitarka A (1998) A technique forsimulating strong ground motion using hybrid Green’s
function. Bull Seism Soc Am 88:357–367
Onishi Y, Horike M (2004) The extended stochastic simulation
method for close-fault earthquake motion prediction andcomments for its application to the hybrid method. J Struct
Constr Eng Trans AIJ 586:37–44
Pitarka A, Somerville P et al (2000) Simulation of near-fault
ground-motion using hybrid Green’s functions. Bull Seism
Soc Am 90:566–586
Tsuruki M, Kagawa T et al (2006) Study on a high-cut filter forstrong ground motion prediction – based on the observed
records during the 1995 Hyogo-ken Nambu Earthquake.
J Jpn Assoc Earthqu Eng 6:94–112 (in Japanese)
Wald DJ, Heaton TH, Hudnut KW (1996) The slip history of
the 1994 Northridge, California, earthquake determinedfrom strong motion, teleseismic, GPS, and leveling data.
Bull Seism Soc Am 86:49–70
Yokoi T, Irikura K (1991) Empirical green's function techniquebased on the scaling law of source spectra. J Seism Soc ofJapan Second Series, 44:109–122 (in Japanese)