Bulletin of the Seismological Society of America, 90, 3, pp. 566–586, June 2000

Simulation of Near-Fault Strong-Ground Motion
Using Hybrid Green’s Functions
by Arben Pitarka, Paul Somerville, Yoshimitsu Fukushima, Tomiichi Uetake,and Kojiro Irikura

Abstract The recently proposed hybrid Green’s function method is designed to
combine the advantages of both deterministic and stochastic approaches to simulating
broadband ground motion when records of small events are not available. The method
has the flexibility of incorporating complexities in the source, wave path, and localsite effects into strong ground motion simulations. In this article we analyze its
effectiveness at simulating near-fault ground motions by comparisons with the empirical source time function method, empirical ground-motion-attenuation relations,and recorded near-fault ground motion. We present a simple model for introducing
the effect of the radiation pattern to the stochastic Green’s functions in the intermediate frequency range (1–3 Hz). The numerical test results of the method and thegenerally good agreement between simulated and recorded ground motion from the
17 January 1995 Kobe earthquake shown in this study indicate that the technique
has the capability of reproducing the main characteristics of near-fault ground motion.
Introduction
broadband strong ground motion simulation if recordings ofsmall earthquakes are not available (e.g., Hartzell et al.,1999). The empirical source time function (ESF) techniqueproposed by Hadley et al. (1982) is one of them. Theirmethod was later improved by Wald et al. (1988) and Somerville et al. (1991) by introducing an empirical representation of the radiation pattern using empirical source functions. The summation of the Green’s functions is done using
Hartzell’s technique (1978). The method proved to be very
effective for simulating ground motions from large earthquakes (e.g., Cohee et al., 1991; Somerville, 1993; Graves
et al., 1996).
Recently Kamae et al. (1998) proposed a new hybrid
technique for simulating near-fault ground motion, which
combines deterministic and stochastic approaches to generate synthetic Green’s functions and uses composite-source
modeling. The stochastic approach used to simulate thehigh-frequency part of the Green’s function is based on theBrune omega-squared point-source model. The strongground motion from a finite fault is then calculated by summing up the synthetic Green’s functions representing the
response of subfaults using Irikura and Kamae’s (1994) summation technique. As demonstrated by Atkinson and Silva(1997), this finite-fault model correctly matches the observed spectral shapes and amplitude. The technique ofKamae et al. (1998) was successfully applied to simulating
ground motions from the Kobe, Japan, earthquake. How-The challenge of near-fault broadband simulation techniques is to provide ground motions with characteristics thatmatch recorded ground motions affected by various source
process and underground structure conditions. Recent analyses of ground motions from damaging earthquakes showthat near-fault ground motions are characterized by longperiod pulses caused by the forward rupture directivity ofthe source and amplified as the waves propagate into adjacent basins, which are often fault controlled. The period,
duration, and amplitude of such pulses has a significant impact on the response of buildings during strong earthquakes.
Most of the long-period characteristics of near-fault motions
can be well simulated using deterministic methods based onkinematic slip models and heterogeneous crustal models. On
the other hand, stochastic models (e.g., Hanks and McGuire
1981; Joyner and Boore, 1986; Silva et al., 1995; Beresnev
and Atkinson, 1998) do not provide a good representationof near-fault ground motions in time domain, but they perform well at predicting the amplitude and frequency content
of motion, on average, at long distances and over a wide
frequency range. Chin and Aki (1991) extended the randomvibration technique to include the effects of the finite-fault
rupture and site-amplification factors. However, phase effects due to the propagation path and site response are notincluded in their approach.
Hybrid techniques that combine both deterministic andempirical approaches can increase the effectiveness of
566

567

Simulation of Near-Fault Strong-Ground Motion Using Hybrid Green’s Functions

ever, their representation of the radiation pattern effect waslimited to the low-frequency range. In this article we introduce a simple technique that extends the effect of the radiation pattern to the intermediate-frequency range (1–3 Hz).
We analyze the effectiveness of the hybrid Green’s function
technique at simulating near-fault ground motions by comparisons with the ESF method, empirical ground-motionattenuation relations, and recorded near-fault ground motion.

Methodology
The hybrid Green’s function (HGF) method is an extension of the semiempirical Green’s function method. It has
been described in detail by Kamae et al. (1998). Here weoutline some important aspects and improvements of thetechnique for generating hybrid Green’s functions. The HGF
method is based on the composite-source-modeling technique. The mainshock fault plane is divided into subfaults,
each of them representing the rupture area of a small event.
The number of subfaults N is determined from equation (1)assuming that the mainshock and the small event have thesame stress drop:N ⳱ (M0m /M0s)2/3(1)

where M0m and M0s are the seismic moments for the mainshock and the small event, respectively.
The mainshock seismogram U(X, t) at a vector position
X is calculated by summing up the synthetic seismogramsu(X, t) representing the response of each subfault as follows(Irikura and Kamae, 1994):N

U(X, t) ⳱

兺

i⳱1

Ci

r
F(t ⳮ ti) * u(X, t)ri

Kamae et al. (1998) tested the technique using fault
models with the stress-drop ratio fixed at a constant value
within the asperities and zero elsewhere on the fault. Herewe consider a more general case where large scale asperities
are introduced by varying the stress drop over the fault surface. The stress drop ratio Ci that appears in equation (2) iscalculated using the approach proposed by Frankel (1995).
Frankel’s technique satisfies the condition that the highfrequency spectral level of the Green’s function sum should
equal that of the mainshock. Also the root mean square stressdrop over the mainshock fault plane should be equal to thestress drop of the small event Drs:

冢

N

Am

1/2

⳱ Drs

(3)

where As and Am are the areas of the subfault and main event,respectively.
Often the source model for a scenario earthquake is represented by the final slip distribution on the fault. To derive
the stress drop-ratio distribution Ci from the slip distribution
di on the fault we assume that they are proportional:
Ci ⳱ ␣ di /dmax(4)where dmax is the maximum slip; ␣ can be determined bysubstituting equation (4) into equation (3) and rearranging:

␣⳱

(2)

冣

兺 Dri2Asi⳱1冢 冢 冣冣
AmAsN
di 2i⳱1 dmax

1/2

兺

(5)

where:Ci ⳱ Dri /Drsti ⳱ (ri ⳮ ro)/b Ⳮ ni /Vr
F(t) ⳱ d(t) Ⳮ1n(Nⳮ1)n兺 d 冤t ⳮ (k ⳮ 1) (N ⳮ 1)n冥k⳱1s

Ci is the ratio of the stress drop Dri in the subfault i and thatof the small event Drs, r, ri, and ro are the distances fromthe site to the hypocenter of the small event, from the site
to ith subfault, and from the site to the starting point ofrupture on the fault, respectively. ni is the distance between
the starting point and the ith subfault. b, Vr, and s are shearwave velocity, rupture velocity, and the rise time of the large
event, respectively. n is an arbitrary integer number used toshift the artificial periodicity to a frequency higher than thatof interest. A stochastic component is included in the speedof the rupture propagation to simulate the heterogeneities in
the kinematics of the rupture.

The procedure we used to calculate hybrid Green’s
functions is similar to the original one in that it computes
the low- and high-frequency ranges separately and then combines the two of them to produce a single time history. The
low-frequency part of the HGF is generated deterministicallyusing numerical techniques that take into consideration the
wave propagation and site effects. The high-frequency partis calculated by a stochastic approach using random vibration time series with a Brune (1970) far-source spectrum anda frequency-dependent radiation-pattern model.

Radiation-Pattern Model
While at low frequencies the radiation pattern of waves
radiated from the seismic source is readily observed; at highfrequencies the pattern is not as clear. Based on analysis ofpeak accelerations from the Whittier Narrows, California,
earthquake and an aftershock, Vidale (1989) argued that forthis earthquake the standard double-couple radiation pattern

568

A. Pitarka, P. Somerville, Y. Fukushima, T. Uetake, and K. Irikura

is observable in the frequency range of 3–6 Hz. On the other
hand, analysis of the ground motion from the 1979 Imperial
Valley earthquake and an aftershock (Liu and Helmberger,
1985) show a double-couple pattern at frequencies smaller
than 1 Hz but not at 2 Hz or higher. The frequency-dependence of the radiation pattern is probably due to scatteringfrom small-scale heterogeneities in the fault zone and aroundthe recording station that tend to affect a broad range of
frequencies. The distribution of these heterogeneities is assumed to be random, and the degree of randomness is assumed to increase with frequency.
In order to account for the frequency-dependent radiation-pattern effect, we have used average radiation coefficients that produce a gradual decay of the radiation patterneffect as a function of frequency (Kamae and Irikura, 1992).
We incorporated the frequency-dependent radiation pattern
into the computation of the HGF. In the tests shown in thisarticle the radiation pattern at frequencies lower than 1 Hz
is modeled using theoretical radiation coefficients of bodywaves for a double-couple point source. At the intermediate
frequency range 1–3 Hz we used average radiation coefficients for body waves. They are calculated using the technique proposed by Boore and Boatwright (1984) assuming
a variable range of takeoff angles that is modeled as a linear
function of frequency. This approach works under the assumption that the waves that carry the energy contributing
to the peak acceleration are body waves coming directlyfrom the source and have takeoff and azimuth angles randomly distributed within given ranges. At frequencies higher
than 3 Hz we used a constant average radiation coefficient
calculated at that frequency. Figure 1 illustrates the frequency-dependent radiation pattern of SH waves for a point
source with pure strike-slip mechanism calculated with theproposed approach.
The choice of the frequency limits 1 Hz and 3 Hz usedin this study is based on the fact that the deterministic approach to modeling underground structure and the source
process is usually reliable at low frequencies only (less than
1 Hz), and on the frequency dependence of the coherence ofthe radiation pattern analyzed by Liu and Helmberger(1985). In the following section we describe the technique
we used to calculate average body-wave radiation coefficients.

Average Body-Wave Radiation Coefficients
Most broadband simulations that apply random vibration theory use the S-wave rms average radiation coefficient.
This average is independent of the fault orientation and ingeneral is calculated over the whole focal sphere, whichmakes the near-fault ground-motion modeling inaccurate at
low frequencies. In our application only a portion of the
focal sphere is sampled using a range of appropriately chosen takeoff angles. The radiation coefficients of SV and SH
waves are averaged over a range of azimuth and takeoff

Figure 1. Frequency-dependent SH-wave radiation pattern for a point source with a pure strike-slip
mechanism and an azimuth angle of 0 degrees. Notethat at frequencies higher than 3 Hz the radiation effect is assumed as purely stochastic; therefore, theradiation pattern transforms into a circle.
angles randomly distributed using Boore and Boatwright’s
technique (1984).
Following their technique the radiation coefficient for agiven type of wave can be calculated using the expression
p 2p

具G典 ⳱

冮 冮 [W(h, ␾)G(h, ␾)] sin (h) d␾dh
0

0

(6)

p 2p

冮 冮 [W(h, ␾)] sin(h)d␾dh
00where ␾ is the azimuth and h is the takeoff angle. It is assumed that the weighting function W(h, ␾) is nonzero andequal to one for takeoff angles h in the range h1 ⬍ h ⬍ h2
and azimuth angles ␾ in the range of (␾r ⳮ ␾1) ⬍ ␾ ⬍ (␾rⳭ ␾1) where ␾r is the receiver azimuth.
The integral in equation (10) is evaluated using a Monte
Carlo scheme (Boore and Boatwright, 1984). According tothat scheme the average value of the radiation coefficient G
can be evaluated as the average coefficient evaluated at M
pairs of hi and ␾i
M

具G典 ⳱

1
G(hi, ␾i)
M i⳱1

兺

(7)

where hi, ␾i are calculated based on M random numbers gi,
vi between 0 and 1 that are mapped as follows:Simulation of Near-Fault Strong-Ground Motion Using Hybrid Green’s Functions

569

Figure 3.

Fault and receiver locations. Stars anddiamond indicate the small events and the mainshock
epicenters, respectively. The triangles indicate the locations of receivers.

Figure 2.

Two different stress-drop-ratio models
used to simulate broadband time histories of theground motion from a M 6.5 earthquake using synthetic Green’s functions from a M 4.6 small event.

hi ⳱ cosⳮ1 [(1 ⳮ ni) cos h1 Ⳮ ni cos h2]
␾i ⳱ ␾r Ⳮ ␾1(0.5 ⳮ gi).

comparisons of synthetic and recorded seismograms, andempirical attenuation relations of strong-ground motion.

(8)

The frequency dependence of the radiation coefficients is
introduced by expressing h1, h2, and ␾1 as linear functionsof frequency between f 1 and f 2:
h1 ⳱ hr ⳮ p/6 ( f ⳮ f1)/( f2 ⳮ f1)
h2 ⳱ hr Ⳮ p/6 ( f ⳮ f1)/( f2 ⳮ f1)
␾1 ⳱ 2p/3 ( f ⳮ f1)/( f2 ⳮ f1)(9)where h1 is the receiver takeoff angle.
The average radiation coefficients at frequencies higher
than f2 are assumed to be equal to those calculated at thefrequency f ⳱ f 2. In the numerical tests of the techniquepresented in this study we used f 1 ⳱ 1 Hz and f 2 ⳱ 3 Hz.

Tests of the Method
In this section we present results of several tests thatwere performed to check the validity of the HGF technique
at simulating near-fault ground motion. The tests consist ofcomparisons with the new ESF method described earlier,Comparison with the ESF Method
We compared synthetic seismograms calculated with
the HGF and ESF methods for a M 6.5 earthquake with apure strike-slip mechanism and a seismic moment M0 ⳱
6.31*1025 dyne cm. The representation of the radiation pattern in the ESF method is done using empirical source functions (Somerville et al., 1991). The method has been tested
against strong-motion data from several types of earthquakes
and has proven to be very effective for simulating broadband
ground motions from large earthquakes (e.g., Cohee et al.,
1991; Somerville et al., 1991). Figure 2 shows the stressdrop models used to simulate broadband ground motions at
stations regularly distributed around the fault. The fault isburied at 1-km depth, and its length and width are 20 km
and 9 km, respectively. Its location and the station configuration are shown in Figure 3. We used a rupture velocityof 2.7 km/sec and a bell-shaped source-time function with aduration of 0.8 sec. In order to avoid additional complexities
in the simulated waveforms due to wave-propagation effects
the simulations were done with an homogeneous half-spacevelocity model with Vp ⳱ 5.2 km/sec, Vs ⳱ 3.0 km/sec, r0
⳱ 2.3 g/cm3 and Q ⳱ 500. The fault plane was divided into
13 ⳯ 5 subfaults each having an area of 1.54 ⳯ 1.8 km2.

570

A. Pitarka, P. Somerville, Y. Fukushima, T. Uetake, and K. Irikura

The number of subfaults was derived using the scaling relationship expressed by equation (1). In the Green’s function
simulations we used six subevents regularly distributed on
the fault plane. Each of them is of magnitude M 4.6 withseismic moment M0 ⳱ 1.26*1023 dyne cm.
The low frequency part (⬍ 1 Hz) of the six hybrid
Green’s functions was calculated using the 1D frequencywave number technique of Saikia (1994). As in the originalscheme of Kamae et al. (1998) the response of each subfault
is represented by the Green’s function corresponding to theclosest subevent corrected for the rupture arrival time andgeometrical spreading.
In Figures 4 and 5 we compare the distribution of thefault-normal peak acceleration and peak velocity, respectively, calculated with the HGF and ESF methods using thesource model shown in Figure 2(a). In both simulations thelow- and high-frequency parts of the synthetic Green’s functions are combined by a pair of matching filters with a transition-frequency range at 0.9–1.1 Hz. The methods agree
relatively well. The peak-velocity pattern, which is similar
in both calculations, reflects the source-directivity effect. Because the rupture starting point is located in the southernpart of the fault, the peak velocity is larger at most of thesites located north of epicenter. The difference in peak velocity between the two techniques is on average less than
15% while the difference in peak acceleration is slightly
higher at stations very close to the central part of the fault.
The near-fault peak-acceleration distribution calculated with
HGF is more uniform than that of peak velocity, which isclearly affected by the source directivity. This result is in a

Figure 4.

Map of peak acceleration (fault-normal component) distribution calculated with the HGF method (left panel) and the ESF method (right panel) for a M 6.5
strike-slip earthquake. Also shown are the locations of the fault-trace (solid line), stations (triangles), small-events epicenter (stars), and mainshock epicenter (hexagon).

Simulation of Near-Fault Strong-Ground Motion Using Hybrid Green’s Functions

Figure 5.

571

Same as Figure 4, but for the peak velocity.

good agreement with the observation that source directivity
does not affect the high frequencies.
Comparisons between synthetic acceleration and velocity waveforms calculated with the two techniques are shownin Figure 6 and Figure 7, respectively. The waveforms andthe time duration of the synthetic seismograms are quite
similar. The long-period pulses seen at receivers locatednorth of the fault, and the increase in the time duration of
the ground motion at receivers located south of the epicenter,are a direct consequence of the source process which is similarly represented in the two techniques. The differences seenin the high-frequency waveforms appearing before the S
wave are due to the fact that the high-frequency part of thehybrid Green’s functions is calculated only for the S wave.
The limited number of hybrid Green’s functions usedin this simulation makes the HGF calculation less accurate
at sites affected by backward directivity. At these sites thephases representing the energy contribution from each sub-fault are distinctively lagged in time since the rupture propagates in the opposite direction, away from the stations.
Therefore even small differences in the synthetic Green’s
functions used by the two techniques can have a significant
effect on the synthetic waveforms. There is no limitation in
the number of hybrid Green’s functions used in the proposed
technique. Increasing the number of hybrid Green’s functions certainly improves the quality of the simulation, but at
the expense of the computation time. The choice of the numerical technique used to calculate the hybrid Green’s functions as well as their number depends on knowledge of theunderground structure and the degree of the source-model
complexity.
In Figures 8 and 9 we show the comparison of the acceleration- and velocity-response spectra, respectively, calculated at receivers 1 to 7, which are affected by source
directivity. The spectra at these stations compare very wellin the period range 0.1 to 10 sec.

572

A. Pitarka, P. Somerville, Y. Fukushima, T. Uetake, and K. Irikura

Figure 6. Time histories of the synthetic acceleration (in cm/sec2) calculated at stations shown in Figure 3 using the HGF method (left panel) and the ESF method (rightpanel) for a M 6.5 strike-slip earthquake. (a) fault-normal component.

We applied the procedure of Abrahamson et al. (1990)
to measure the goodness of fit of response-spectral acceleration between the ground motions simulated with the two
methods. In our experiment the goodness of fit is measured
by the bias, which is expressed as a function of period. The
bias measures the difference between the two simulations

averaged over all stations at the considered period range.
Figure 10 shows the bias and the standard error of residual
between the response spectra simulated with the two techniques. The goodness of fit in response-spectral accelerationof the two horizontal components indicates that the bias
mean between the two simulations is relatively small, al-Simulation of Near-Fault Strong-Ground Motion Using Hybrid Green’s Functions

Figure 6.

573(Continued) (b) fault-parallel component.

though there is a tendency for the HGF method to produce
slightly larger ground motions. The bias at long periods andthe scattering are greater at stations that see the whole fault
plane and are not affected by the rupture directivity. As discussed previously the bias at these stations is attributable tothe difference in the number of Green’s functions used by
the two techniques.

Comparisons with Empirical Response Spectral
Acceleration Models
The performance of the HGF technique was tested bycomparing simulated- and empirical-response spectra forrock sites using the Abrahamson and Silva (1997) model for
strike-slip earthquakes of M 6.5. The empirical relation of

574

A. Pitarka, P. Somerville, Y. Fukushima, T. Uetake, and K. Irikura

Figure 6.

(Continued) (c) vertical component.

Abrahamson and Silva for crustal earthquakes was derived
using a worldwide data set of shallow events. According totheir soil classification, sites with shear-wave velocity higher
than 600 m/sec are considered as rock. We chose M 6.5 forthe comparison because it is constrained by the largest
amount of data.
Figure 11 shows the comparisons of simulated and empirical peak-ground-acceleration attenuation, and spectral-

acceleration attenuation at several periods. The synthetic accelerations were obtained based on the simulation describedin the previous section (white circles) using the source model
shown in Figure 2(a) and a second simulation (black circles)
using the source model shown in Figure 2(b). In general thesimulated accelerations fall within the range of one standard
deviation in the empirical model for all the rupture distancesand periods except for 2 sec. At this period the simulations

Simulation of Near-Fault Strong-Ground Motion Using Hybrid Green’s Functions

Figure 7.

575

Same as Figure 6, but for the velocity (units are cm/sec).

give slightly larger amplitudes at some of the stations locatedwithin the rupture distance range of 5–20 km. We attribute
this discrepancy to the directivity effect, which is particularly strong at stations located north of the epicenter in thecase of the first simulation and above the hypocenter in thecase of the second simulation. The effect of the directivity
on simulated ground motions is also seen when we compare
the results from the two source models (see Fig. 11). Bothmodels give comparable peak accelerations at periods
shorter than 1 sec. At periods of 1 and 2 sec the first model
produces stronger directivity giving slightly larger amplitudes at most of the stations.
The same source representation and matching filters, described in the previous section, were used to compare thesimulated average horizontal spectral accelerations and the
median predicted by Abrahamson and Silva’s (1997) model

576

A. Pitarka, P. Somerville, Y. Fukushima, T. Uetake, and K. Irikura

Figure 7.

(Continued).

for rock sites at a distance from the rupture of 10 km (Fig.
12). For this comparison, we generated synthetic ground motion at 11 stations located around the fault with a distancefrom the rupture of approximately 10 km. We used the 1D
velocity model shown in Table 1 to calculate the synthetic
Green’s functions. In this model the shear-wave velocity inthe upper 100 m is 1.2 km/sec. The synthetic ground motions
were calculated with three different source models, allowing

us to investigate the effect of the source on the simulated
spectral acceleration at long periods. Although the empiricalresponse spectral acceleration is based on data from a wide
range of strike-slip earthquakes, the comparison with our
synthetic spectra derived from a single realization for eachsource model is good, especially at periods smaller than 1
sec. This is partly due to the fact that the high-frequency partof the near-fault ground motion is less sensitive to the large-

Simulation of Near-Fault Strong-Ground Motion Using Hybrid Green’s Functions

Figure 7.

577(Continued).

scale complexities in the rupture process. The peak at around
0.2 sec and the quick fall-off of the spectrum at periods
shorter or longer than this value are matched very well. At
periods longer than 1 sec the simulations give higher response-spectral accelerations. The comparison between thethree models suggest that the response at these periods ismainly controlled by the large-scale features of the source
model. Model II and Model III that contain multiple

smoothed asperities give better results. The relatively large
and smooth asperity in Model I, the location of the rupture
starting point, and the pure strike-slip mechanism enhancethe directivity effect, which is strong at periods longer than
1 sec at most of the stations. Numerical experiments (e.g.,Graves, 1998) performed with several source models suggestthat rough and multiple asperities tend to diminish the coherency of long-period energy release during the fault rup-

578

A. Pitarka, P. Somerville, Y. Fukushima, T. Uetake, and K. Irikura

Figure 8.

Comparison of acceleration response spectra calculated at seven stations
using the HGF method (solid line) and the ESF method (dashed line).

Figure 9.

Same as Figure 8, but for the velocity response spectra (5% damping).

Simulation of Near-Fault Strong-Ground Motion Using Hybrid Green’s Functions

579

Figure 10. Upper panels: Model bias (solid line) of the natural logarithm of spectral
acceleration (reflecting average ratio of synthetic ground motion simulated with HGF
and ESF methods), and its 90% confidence interval (dashed lines) at stations shown in
Figure 3. A positive bias indicates that the HGF method gives larger ground motion
than the ESF method. Lower panel: Standard error of residual between the response
spectra obtained with the two techniques, with and without correction for model bias,
for the fault-normal component (left panels) and fault-parallel component (right panels).

ture. Consequently they produce weaker near-fault groundmotions.
Our result suggests that in order to increase the effectiveness of the HGF technique at simulating ground motions
from scenario earthquakes it is necessary to systematicallyspan the range of possible slip models. The departure fromthe empirical model at stations affected by rupture directivity
can be significant for smooth asperity source models.

Simulation of Near-Fault Ground Motion from the1995 Kobe Earthquake
We applied the HGF method to simulate records of nearfault ground motion from the 1995 Kobe earthquake. The
location of the fault and five stations considered in this study

are shown in Figure 13. Among these stations, only KBU isa rock site (Table 2). The basin depth varies from about 600
m beneath KOB and MOT, which are located very close to
the basin edge, to 1.3 km beneath TKT and FKA. A schematic view of the fault model used in the simulation is shownin Figure 14. The model has the following features, proposedby Kamae et al. (1998). It includes four square-shaped asperities. The stress drop is assumed to be variable withineach asperity and zero in the other regions of the fault. The
location of the asperities roughly corresponds to the areas oflarge slip concentration derived from inversions of strongmotion data (e.g. Sekiguchi et al., 1996; Wald, 1996). We
found that small modifications of the fourth asperity area do
not affect significantly the ground motion at sites consideredhere. This can be explained by the fact that the sites are580

A. Pitarka, P. Somerville, Y. Fukushima, T. Uetake, and K. Irikura

Figure 11. Comparisons of simulated (circles) and empirical peak ground acceleration attenuation and spectral acceleration attenuation (solid line) of Abrahamson andSilva (1997), for a M 6.5 strike-slip earthquake. Dashed line is the standard deviationof the empirical attenuation. The simulated accelerations were obtained using the source
models shown in Figure 2(a) (white circles) and Figure 2(b) (black circles).

located outside the area affected by the forward-source directivity caused by the rupture of the fourth asperity.
Hybrid Green’s functions are calculated for four point
sources located at the center of each asperity assuming amagnitude 4.7 for each source. The source parameters for
the M 4.7 small event are summarized in Table 3. We used
1D and 3D techniques to simulate the low-frequency part
(0.1–1 Hz) of the Green’s function and the stochastic technique of Boore (1983) for the high-frequency part (⬎1 Hz).
The transition-frequency range in the matching filters was
0.9–1.1 Hz.
Pitarka et al. (1998) found that, at long periods, the local
bedrock geometry does not have a significant effect at KOB
and MOT, while it has an impact at TKT and FKA. According to Kamae et al. (1998), who tested the 3D versus
1D effect at the basin sites considered in this study, the effectof the thin shallow sedimentary layers is significant at frequencies around 1 Hz. Based on these results, we used 1D

velocity models to calculate the low-frequency part of theGreen’s functions at KBU, KOB, and MOT (Table 4), andthe 3D velocity model of Pitarka et al. (1998) for TKT and
FKA. In the 3D calculations we used a variable-grid finitedifference scheme (Pitarka, 1999). The velocity model was
discretized with a grid spacing of 60 m in the horizontal
direction, and 60 m, 100 m, 200 m, 300 m, and 400 m inthe depth intervals of 0–1 km, 1–2 km, 2–5 km, 5–8 km,
and 8–24 km, respectively.
The simulated and observed velocity and acceleration
seismograms at KBU, KOB, and MOT are depicted in Figure 15. Except for the fault-normal component of acceleration at KOB, there is a good agreement in the peak velocity
and acceleration at all sites. The agreement between the simulated and recorded acceleration waveforms is fair and thatof the velocity waveforms is relatively good. The observed
long-period pulses are well simulated.
The comparison at TKT and FKA is shown in Figure

581

Simulation of Near-Fault Strong-Ground Motion Using Hybrid Green’s Functions

Figure 12. Comparisons of simulated (thin line) and empirical (thick line) spectral
acceleration of Abrahamson and Silva (1997) at a rupture distance of 10 km for a M
6.5 strike-slip earthquake. The simulations were done using three different source models shown in the lower panels.

Table 1
Velocity Model for Rock Sites
Depth (km)Vp (km/sec)Vs (km/sec)

Density (g/cm3)

Qp

Qs

0.10
1.25
4.45
⬎4.45

2.4
3.5
4.4
5.7

1.2
2.0
2.5
3.5

2.0
2.1
2.2
2.3

100
200
500
500

50
100
300
300

16. It is obvious that the complexities in the waveforms
observed at these sites are related not only to the source
process, but also to the basin-edge effect and the surfacesedimentary-layer response. Our simulation reproduces well
the duration and the first part of the ground motion, which
is mainly affected by the rupture directivity and to a lesser
extent by the basin edge effect. We interpret the second large
pulse seen in both sites as a superposition of waves coming

directly from the fourth asperity and surface waves generated at the basin edge. The later phases in the acceleration
observed at TKT suggest that the response of the shallow
sediments during the mainshock was nonlinear. They were
not reproduced by our elastic technique. Aguirre and Irikura(1997), who modeled the response of reclaimed land sites in
the Kobe region, concluded that the response of soft sediments throughout the damaged region may have been nonlinear during the mainshock.

Conclusions
The HGF method is designed to combine the advantagesof both deterministic and stochastic approaches and has theflexibility of incorporating complexities in the source, wavepath, and local site effects into strong ground motion simulation. It is sometimes argued that the effects of the sourceand path must be averaged over various scenarios in pre-

582

A. Pitarka, P. Somerville, Y. Fukushima, T. Uetake, and K. Irikura

Figure 13. Map of Kobe area and station locations. The rectangle shows the areal
extent of the basin model used in the 3D simulations of the Green’s functions. Also
shown are the locations of the damaged zone and epicenters of four hypothetical small
events (stars).

Table 2
Station Locations
Station Name

Lat. (deg.)Long. (deg.)

Geological Conditions

KBU
KOB
MOT
TKT
FKA

34.72534.68834.72534.64934.695135.240135.180135.281135.138135.211

Soft Rock
Alluvium
Alluvium
Alluvium
Alluvium

dicting ground motion for seismic design, so the HGF
method may not perform better than simpler techniques that
represent the time history of the ground motion using thestochastic approach. This is not quite true in the case of nearfault ground motion, especially that affected by source directivity and basin-edge geometry. Records of near-fault
ground motion from 1989 Loma Prieta, 1992 Landers, 1994
Northridge, and 1995 Kobe earthquakes demonstrate thatmost of the seismic energy is carried by long-period waves
generated by source directivity and affected by the wave
path (e.g., Steidl et al., 1991; Wald et al., 1991; Cohe andBeroza, 1994; Olsen and Archuleta, 1996; Somerville et al.
1996; Kamae and Irikura, 1998). Such ground motions cannot be well represented by stochastic models in time domain.
The long-period characteristics of near-fault ground motion
have a significant impact on the dynamic response of large
structures. Current research on the dynamic response ofstructures to near-fault ground motions relies heavily on
nonlinear dynamic analysis using ground-motion time histories as input. This requires the use of several strong-motion
seismograms, recorded or simulated, each of which is therealization of a single scenario that incorporates particularsource and wave-path effects at low frequencies. For thispurpose it is not meaningful to use a ground-motion time
history that represents the average of multiple scenarios. The
HGF method is designed to fulfill such a requirement, including the effect of the coherence of the source radiation
pattern, which is at the origin of the source-directivity effect.
We believe that the period, amplitude, and timing of thelong-period phases that often characterize near-fault groundmotion are of primary importance for building design, andthat the improvement of numerical techniques and knowledge of underground structure are needed for the development of methods that apply a deterministic approach to nearfault ground-motion prediction.
Although it occupies a narrow-frequency range, the inclusion of the frequency-dependent radiation-pattern effectinto broadband near-fault ground-motion simulations provides an important enhancement of the method. The simplistic representation of the radiation-pattern effect by constant coefficients usually applied in broadband simulations
may produce unrealistic synthetic seismograms. Studies of
the radiation-pattern effect in a broad-frequency range basedon ground-motion data are necessary to understand the frequency dependence of the source radiation and its effect innear-fault ground motion.

583

Simulation of Near-Fault Strong-Ground Motion Using Hybrid Green’s Functions

Figure 14. Source model consisting of four rectangle-shaped asperities. The stressdrop ratio within the asperities is estimated by the forward modeling of near-fault
ground motion.

Table 3
Source Parameters for M 4.7 Small Event
Seismic Movement (dyne cm)
Corner Frequency (Hz)
Focal Mechanism
Rise Time (sec)

1.30 ⳯ 1023
1.3
strike 53⬚; dip 90⬚; rake 180⬚
0.12

Table 4
One-Dimensional Velocity Models
Site

Depth
(km)

Vp
(km/sec)

Vs(km/sec)

Density(g/cm3)

Qp

Qs

KBU
KBU
KBU
MOT
MOT
MOT
MOT
MOT
MOT
MOT
KOB
KOB
KOB
KOB
KOB
KOB

0.5
5.0
⬎5.0
0.035
0.05
0.15
0.30
0.6
5.0
⬎5.0
0.06
0.2
0.4
0.6
5.0
⬎5.0

3.2
5.7
6.0
1.5
1.7
1.8
1.8
1.9
5.7
6.0
1.5
1.7
1.8
1.9
5.7
6.0

1.8
3.2
3.46
0.3
0.45
0.50
0.65
1.1
3.2
3.46
0.45
0.55
0.65
1.1
3.2
3.46

2.0
2.4
2.7
1.9
1.9
2.0
2.0
2.2
2.4
2.7
1.9
2.0
2.0
2.2
2.4
2.7

500
1000
1000
130
150
150
150
300
1000
1000
150
150
150
300
1000
1000

300
400
500
50
80
80
80
170
400
500
80
80
80
170
400
500

The generally good agreement between simulated andrecorded ground motion obtained for the 1995 Kobe earthquake and shown in this study indicates that the technique
has the capability of reproducing the main characteristics ofnear-fault ground motion from basin-edge faults. At sites
where the nonlinear effect of sediments is preponderant, thetechnique needs to be augmented by nonlinear techniques.
Numerical experiments indicate that the summation ofGreen’s functions using the technique proposed by Irikura
and Kamae (1994) performs relatively well when the seismic
source is characterized by well-distinguished large asperities, as it is the case of the Kobe, Japan, earthquake.

Acknowledgments
Discussions with K. Kamae and T. Sato are greatly appreciated. We
thank G. Atkinson and an anonymous reviewer for their valuable suggestions. This work was partially supported by TEPCO, Japan.

References
Abrahamson, N. A., P. G. Somerville, and C. Allin Cornell (1990). Uncertainty in numerical strong motion predictions, in Proc. 4th Nat. U.S.
Conference Earthquake Engineering (Palm Springs), Earthquake
Eng. Res. Inst. 1, 407–416.
Abrahamson, N. A., and W. J. Silva (1997). Empirical response spectral
attenuation relations for shallow crustal earthquakes, Seism. Res. Lett.
68, 94–127.
Aguirre, J., and K. Irikura (1997). Nonlinearity, liquefaction, and velocity

584

A. Pitarka, P. Somerville, Y. Fukushima, T. Uetake, and K. Irikura

Figure 15. Comparison of synthetic (upper trace) and recorded (lower trace) time
histories of the ground motion at KBU, KOB, and MOT. (a) Velocity: the number tothe right of each trace gives its maximum amplitude in cm/sec (b) acceleration: thenumber to the right of each trace gives its maximum amplitude in cm/sec2.

Simulation of Near-Fault Strong-Ground Motion Using Hybrid Green’s Functions

Figure 16.

585

Same as Figure 15, but for TKT and FKA.

variation of soft soil layers in Port Island, Kobe, during the Hyogoken Nanbu earthquake, Bull. Seism. Soc. Am. 87, 1244–1258.
Atkinson, G. M., and W. Silva (1997). An empirical study of earthquake
source spectra for California earthquakes, Bull. Seism. Soc. Am. 87,
97–111.
Beresnev, I. A., and G. M. Atkinson (1998). Stochastic finite-fault modelingof ground motions from the 1994 Northridge, California, earthquake.
I. Validation on rock sites, Bull. Seism. Soc. Am. 88, 1392–1401.
Boore, D. M. (1983). Stochastic simulation of high-frequency ground motions based on seismological models of the radiated spectra, Bull.
Seism. Soc. Am. 73, 1865–1894.
Boore, D. M., and J. Boatwright (1984). Average body-wave radiation coefficients, Bull. Seism. Soc. Am. 74, 1615–1621.
Brune, J. (1970). Tectonic stress and the spectra of seismic shear waves
from earthquakes, J. Geophys. Res. 75, 4997–5009.
Chin, B. H., and K. Aki (1991). Simultaneous study of the source, path andsite effects on strong ground motion during the 1989 Loma Prieta
earthquake: a preliminary result on pervasive nonlinear site effects,
Bull. Seism. Soc. Am. 81, 1859–1884.
Cohee, B. P., and G. C. Beroza (1994). Slip distribution of the 1992 Landers

earthquake and its implications for earthquake source mechanics,
Bull. Seism. Soc. Am. 84, 692–712.
Cohee, B. P., P. G. Somerville, and N. A. Abrahamson (1991). Simulated
ground motion for hypothesized Mw ⳱ 8 subduction earthquakes in
Washington and Oregon. Bull. Seism. Soc. Am. 81, 28–56.
Frankel, A. (1995). Simulating strong motion of large earthquakes using
recordings of small earthquakes: the Loma Prieta Mainshock as a test
case, Bull. Seism. Soc. Am. 85, 1144–1160.
Graves, R. W. (1998) 3D finite-difference modeling of the San Andreas
Fault: source parameterization and ground motion levels. Bull. Seism.
Soc. Am. 88, 881–897.
Graves, R. W., P. G. Somerville, and N. F. Smith (1996). Simulated broadband strong motion time histories throughout Osaka Prefecture forfive faults beneath and bounding Osaka basin, including near-fault
rupture directivity and 3D basin response effects, Report to Georesearch Institute, Osaka, Japan.
Hadley, D. M., D. V. Helmberger, and J. A. Orcutt (1982). Peak acceleration scaling studies, Bull. Seism. Soc. Am. 72, 959–979.
Hanks, T., and R. McGuire (1981). The character of high-frequency strong
ground motion, Bull. Seism. Soc. Am. 71, 2071–2095.

586

A. Pitarka, P. Somerville, Y. Fukushima, T. Uetake, and K. Irikura

Hartzell, S. H. (1978). Earthquake aftershocks as Green’s functions, Geophys. Res. Lett. 5, 1–4.
Harztell, S., S. Harmsen, A. Frankel, and S. Larsen (1999). Calculation ofbroadband time histories of ground motion: comparison of methods
and validation using strong-ground motion from the 1994 Northridge
earthquake, Bull. Seism. Soc. Am. 89, 1484–1504.
Irikura, K., and K. Kamae (1994). Estimation of strong ground motion inbroad-frequency band based on a seismic source scaling model andan empirical Green’s function technique, Annali di Geofisica xxxvii,
6, 1721–1743.
Joyner, W. B., and D. M. Boore (1986). On simulating large earthquakes
by Green’s function addition of smaller earthquakes, in Earthquake
Source Mechanics, Maurice Ewing Series 6, S. Das, J. Boatwright,
and C. H. Scholz (Editors), American Geophysical Union, Washington, D.C., 37, 269–274.
Kamae, K., and K. Irikura (1992). Prediction of site-specific strong groundmotion using semiempirical methods, in Earthquake Engineering,
Tenth World Conference, 1992, Madrid, Spain.
Kamae, K., and K. Irikura (1998). Source model of the 1995 Hyogo-ken
Nanbu earthquake and simulation of near-source ground motion, Bull.
Seism. Soc. Am. 88, 400–412.
Kamae, K., K. Irikura, and A. Pitarka (1998). A technique for simulatingstrong ground motion using Hybrid Green’s function, Bull. Seism.
Soc. Am. 88, 357–367.
Liu, H. L., and D. V. Helmberger (1985). The 23:19 aftershock of the 15
October 1979 Imperial valley earthquake: more evidence for an asperity, Bull. Seism. Soc. Am. 75, 689–708.
Olsen, K. B., and R. J. Archuleta (1996). Three-dimensional simulation ofearthquakes on the Los Angeles faults system, Bull. Seism. Soc. Am.
86, 575–596.
Pitarka, A. (1999). 3D elastic finite-difference modelling of seismic motion
using staggered-grid with non-uniform spacing. Bull. Seism. Soc. Am.
89, 54–68.
Pitarka, A., K. Irikura, T. Iwata, and H. Sekiguchi (1998). Three-dimensional simulation of the near-fault ground motion for the 1995 Hyogoken Nanbu (Kobe), Japan, earthquake, Bull. Seism. Soc. Am. 88, 428–
440.
Saikia, C. K. (1994). Modified frequency-wavenumber algorithm for regional seismograms using Filon’s quadrature: modelling of Lg waves
in eastern North America, Geophys. J. Int. 118, 142–158.
Sekiguchi, H., K. Irikura, T. Iwata, Y. Kakehi, and M. Hoshiba (1996).
Minute locating of fault planes and source process of the 1995 Hyogoken-Nanbu earthquake from the waveform inversion of strong
motion. J. Phys. Earth 44, 473–487.
Silva, W. J., I. G. Wong, and R. B. Darragh (1995). Engineering characterization of earthquake strong ground motions in the Pacific Northwest, in Assessing Earthquake Hazards and Reducing Risk in the
Pacific Northwest, U.S. Geol. Surv. Profess. Pap. 1560.
Somerville, P. G. (1993). Engineering applications of strong ground motion
simulation, Tectonophysics 218, 195–219.
Somerville, P. G., M. Sen, and B. Cohee (1991). Simulation of strong

ground motion recorded during the 1985 Michoacan Mexico and Valparaiso, Chile earthquakes. Bull. Seism. Soc. Am. 81, 1–28.
Somerville, P. G., C. Saikia, D. Wald, and R. Graves (1996). Implicationsof the Northridge earthquake for strong ground motions from thrustfaults, Bull. Seism. Soc. Am. 86, S115–S125.
Steidl, J. H., R. Archuleta, and S. H. Hartzell (1991). Rupture of the 1989
Loma Prieta, California, earthquake. Bull. Seism. Soc. Am. 81, 1573–
1602.
Vidale, J. H. (1989). Influence of focal mechanism on peak accelerationsof strong motions of Whittier Narrows, California, earthquake and anaftershock, J. Geophys. Res. 94, 9607–9613.
Wald, D. J. (1996). Slip history of the 1995 Kobe, Japan earthquake determined from strong motion, teleseismic and geodetic data, J. Phys.
Earth. 44, 489–503.
Wald, D. J., L. J. Burdick, and P. G. Somerville (1988). Simulation ofacceleration time histories close to large earthquakes, in Earthquake
Engineering and Soil Dynamics II—Recent Advances in Ground Motion Evaluation, Geotechnical Special Publication 20, J. Lawrence
Von Thunn, (Editor), American Society of Civil Engineering, New
York, 430–444.
Wald, D. J., D. V. Helmberger, and T. H. Heaton (1991). Rupture modelof the 1989 Loma Prieta earthquake from the inversion of the strongmotion and broad-band teleseismic data, Bull. Seism. Soc. Am. 81,
1540–1572.
URS Greiner Woodward Clyde
566 El Dorado Street, Suite 100
Pasadena, California 91101(A. P., P. S.)
Izumi Research InstituteShimizu Corporation
Fukoku Seimei Building
2-2-2 Uchisaiwai-cho, Chyoda-ku, Tokyo, 100-0011
Japan(Y. F.)
Power Engineering R&D Center,
Tokyo Electric Power Company
4-1, Egasaki-cho, Tsurumi-ku, Yokohama, 230-8510
Japan(T. U.)
Disaster Prevention Research InstituteKyoto University
Gokasho, Uji, 611
Japan(K. I.)
Manuscript received 29 July 1999.