Bulletin of the Seismological Society of America, Vol. 93, No. 3, pp. 1154–1168, June 2003

A Theoretical Method for Computing Near-Fault Ground Motions in
Layered Half-Spaces Considering Static Offset Due to Surface Faulting,with a Physical Interpretation of Fling Step and Rupture Directivity
by Yoshiaki Hisada and Jacobo Bielak

Abstract An efficient mathematical method is presented for computing the nearfault strong ground motions in a layered half-space, giving explicit consideration to
the static offset due to surface faulting. In addition, the combined effects of “flingstep” and “rupture directivity” on the near-fault ground motions are investigated.
First, after checking the fault integration in the representation theorem, it is foundthat when an observation point is close to the fault plane, Green’s functions exhibit
near singularities, which consist of extremely sharp peaks in a narrow band close to
the observation point. Therefore, direct numerical integration becomes quite onerous
for computing near-fault ground motions, because the dynamic Green’s functions
must then be distributed very densely in order to evaluate accurately the effects of
the near singularities. Instead, a new form of the representation theorem is introduced,
which exploits the property that the dynamic Green’s functions can be approximatedby the corresponding static Green’s functions in the vicinity of the singularities. The
modified theorem, which involves the device of adding and subtracting the static
Green’s functions from the dynamic ones, is the sum of two fault integrations. The
first integration involves the difference of the dynamic and the corresponding static
Green’s functions, while the second contains only the static Green’s functions. This
formulation requires much less CPU time than the original one when near-fault
ground motions are considered, because the near singularities of the dynamic Green’s
functions in the first integration are completely eliminated by subtracting the static
Green’s functions. While the second integration does require a densely distributed
set of points to capture the near-singular behavior of the static Green’s function, it
needs to be performed only once, as it is valid for all frequencies. Subtraction of thestatic Green’s functions from the dynamic functions has the added benefit of making
the integration over the wavenumber in the determination of the Green’s functions
much more efficient, especially for surface faulting. This is because the differenceof the dynamic and static integrands converges rapidly to zero with increasing wavenumbers, whereas the original integrands diverge in the case of a source point on thefree surface.
The proposed methodology is used to investigate the two most important effects
in near-fault ground motions, fling step (e.g., Abrahamson, 2001) and rupture directivity (e.g., Somerville et al., 1997), by paying special attention to the contributionof static and dynamic Green’s functions. It is found that the fling effects stem mainly
from the second integral in the modified representation theorem, which involves thestatic Green’s function. The fling effects are dominant in the slip direction only inthe vicinity of the surface fault and are negligible for buried faults, because the static
Green’s function attenuates rapidly with distance from the fault, r, as the order of(1/r2). Also, more importantly, when an observation point is located above a buried
fault, the medium has to remain continuous, and thus cannot “fling.” By contrast, thedirectivity effects stem mainly from the first integral, which involves the dynamic
Green’s function, and attenuate much more slowly than the fling, on the order from1/r to 1/冪r. The directivity effects are dominant in the fault-normal direction, especially
1154

1155

A Theoretical Method for Computing Near-Fault Ground Motions in Layered Half-Spacesin the forward rupture direction, not only for the surface fault, but also for the buried
fault. Due to the combined effects of fling and directivity in the vicinity of the surface
fault, the directions of the maximum velocities and displacements are inclined withrespect to the fault plane. On the other hand, when softer surface layers are added tothe medium, the directivity effects become more significant than the fling effects,
because the dynamic Green’s functions are more pronounced than the static ones.

Introduction
Recent large-scale inland earthquakes, especially the1999 Chi-Chi earthquake of Taiwan and the 1999 Kocaeli
earthquake of Turkey, dramatically demonstrated the strong
effects of surface faulting on near-fault ground motions (e.g.,
Shin and Teng, 2001). The effects of surface faulting due totectonic deformations have been recently called “fling step”(e.g., Abrahamson, 2001). In contrast to “rupture directivity”
effects, which show the long-period pulses in the direction
normal to the fault plane (e.g., Aki, 1968; Archuleta andHartzell, 1981; Somerville et al., 1997), the fling effects exhibit long-period pulses and permanent static offsets in the
direction of the fault slip. In particular, when the fault rupture reaches the free surface, the faulting slip and the subsequent ground motions are strongly amplified (e.g., Archuleta and Frazier, 1978). Thus, the combination of fling anddirectivity can have destructive effects on various lifeline
facilities and long-period structures, such as high-rise orbase-isolated buildings, in the vicinity of the surface fault.
As an example, Figure 1a,b shows the location of thesurface faults of the 1992 Landers earthquake and the strong
motion at a station in the Lucerne Valley. The station is
located in the vicinity of the surface fault (about 2 km away;e.g., Chen, 1995) and in the forward direction of the fault
rupture. Figure 1c,d shows the fault-normal and -parallel
components of the observed velocity and displacement, respectively. The records clearly show the “forward directivity
pulse” in the fault-normal component (i.e., the long-period
pulse in velocity and the smaller static offset in displacement) and also the fling step in the parallel component (i.e.,
the long-period pulse in velocity and the large static offset
in displacement). Interestingly, because of the combined effects of directivity and fling, the directions of the maximum
velocity and displacement are inclined to the surface fault at
about 45⬚, as shown in Figure 1b.
In this article, we propose an efficient method for simulating these near-fault ground motions using an extended
kinematic source. Even though various numerical methods
have been available to reproduce the near-fault ground motions considering the surface faulting (e.g., Archuleta andFrazier, 1978; Oglesby and Day, 2001; Ramancharla and
Meguro, 2001), analytical methods based on the representation theorem remain a useful and powerful tool. As far as
the analytical Green’s functions are concerned, earlier studies used to assume simple media, such as the homogeneousfull-space (e.g., Aki, 1968; Haskell, 1969) or the homoge-

neous half-space (Kawasaki et al., 1973; Kawasaki, 1975;
Anderson, 1976; Levy and Mal, 1976). Nowadays, one canreadily obtain those of flat-layered half-spaces (e.g., Bouchon, 1979; Apsel and Luco, 1983; Luco and Apsel, 1983;
Hisada, 1993, 1995). However, in order to compute the nearfault ground motions in a layered half-space considering surface faulting accurately, we need to overcome two main obstacles related to the singularity of Green’s functions.
We shall first describe these obstacles using the analytical method. The displacement due to a kinematic fault
model can be expressed as follows in the frequency domain:
Uk(Y; x) ⳱

冮 T (X, Y; x)D (X; x)dR,
R

ik

i

(1)

where Uk is the kth component of displacement in the Cartesian coordinate system at an observation point, Y, x is thecircular frequency, X is a source point on the fault plane,
and R is the fault plane. The Tik is the traction Green’s function, and Di is the ith component of the fault slip; we use
the summation convention for subscript i. The Green’s function of the layered half-spaces is obtained using the following wavenumber integration:
Tik(X, Y; x) ⳱⬁冮 t (X, Y; x, k)dk,
0

ik

(2)

where tik is the integrand of the traction Green’s function.
When the observation point is located on the fault plane, theGreen’s functions become singular as the order of (1/r2).
However, in that case, the fault integration in equation (1)can be evaluated analytically considering Cauchy’s principal
value if the material is homogeneous, and its value is half ofthe slip (e.g., Kane, 1994). On the other hand, when the observation point is not on the fault plane but very close to it,
Green’s functions become nearly singular, that is, their amplitudes show sharp peaks in a narrow area close to the observation point. Therefore, the first obstacle is the need to
distribute Green’s functions very densely on the area close to
the observation point in order to evaluate the fault integration
numerically. This means that a large amount of CPU time isrequired when the dynamic Green’s functions of layered halfspaces are used (e.g., Kawasaki, 1975). In addition, when we
assume a surface fault, we need to calculate numerous

1156

Y. Hisada and J. Bielak

Figure 1. (a) Map of California with the site location; (b) the surface faults and theepicenter of the 1992 Landers earthquake, together with the location of the observation
station at the Lucerne valley; (c) the velocities; and (d) the displacements at the station.
Panel (b) also shows the direction of the strike slip, the directions of the fault-normal and-parallel components, and the directions of the maximum velocity and displacement.

Green’s functions with shallow source points. Therefore, thesecond obstacle is that the integrands of wavenumber integrations (equation 2) do not converge with wavenumber
when the depths of source points are close to or on the freesurface (e.g., Apsel and Luco, 1983; Hisada, 1993, 1995).
In particular, the convergence is extremely slow in the caseof the static Green’s function (x ⳱ 0). Therefore, special
techniques are needed to overcome the two obstacles.
The purpose of this article is to propose a mathematicalmethodology for computing near-fault ground motions effectively and to use it for investigating the effects of flingand directivity in several simple situations. We first carefully
check the fault integration (equation 1) using the simplest
fault model: an axially symmetric circular fault in a homogeneous full-space. Based on the results from this simple
case, we will then propose a new form of the representation
theorem for calculating the fault integration efficiently formore general cases, involving arbitrary kinematic faulting
models in layered half-spaces. In addition, we propose an
efficient method for calculating the wavenumber integration(equation 2), considering the surface faulting. Finally, we
check the validity of the proposed method and investigate
the physical basis of the fling and directivity effects.

Efficient Methods for Computing Near-Fault Ground
Motions in Layered Half-Spaces
Near-Fault Ground Motions Using an Axially
Symmetric Fault Model in a Homogeneous
Full-SpaceWe first check the basic characteristics of the dynamic
and static Green’s functions in the fault integration (i.e.,equation 1) to find efficient ways for computing the nearfault ground motions. In this section, we use the simplest
fault model, that is, the axially symmetric circular fault
model in a homogeneous full-space. In addition, we will
check the attenuation relation of the static offset using thismodel.
Figure 2 shows the fault model and the location of anobservation point. R is the radius of the circular fault model.
We assume a uniform slip, D, over the fault plane. The observation point is located at a distance, z, above the centerof the fault. The dynamic displacement, U, in the same direction as D, is easily obtained by substituting Green’s function of the homogeneous full-space (e.g., Kane, 1994) into
equation (1),

1157

A Theoretical Method for Computing Near-Fault Ground Motions in Layered Half-SpaceszD
4

U(z, x) ⳱冮 f 冦2u Ⳮ u 冢f冣 冧dr,
R

r 2

r

0

1

3

(3)

2

where D is the fault slip, andf ⳱ 冪r2 Ⳮ z2 ,(4)

u1 ⳱ fv ⳮ r2w,f, u2 ⳱ 2(f2v,f ⳮ2fv),

(5)

v ⳱

Figure 2. An axial-symmetric fault model in homogeneous full-space. An observation point is located above the center of the circular fault.

3
3
es
p 2 3
3
ep
ⳮ
Ⳮ
1
ⳮ
Ⳮ
Ⳮ
1
,
s2
s
f
s
p2
p
f

冦

冧

冦

v,f ⳱ ⳮ

冢冣 冦

冧

9
9
4
s s
p 2
9
9
4
p
Ⳮ
ⳮ
Ⳮ
1
e
ⳮ
ⳮ 3 Ⳮ 2 ⳮ Ⳮ 1 2 ep,
s3
s2
s
f2
s
p
p
p
f

冧

冢冣 冦
p 2

冧

(6)

冦 s Ⳮ s ⳮ s Ⳮ 1冧 f e ⳮ冢s冣 冦ⳮ p Ⳮ p Ⳮ p冧 f e ,w,f ⳱ ⳮ

3

3

3

2

2

s

s

2

3

3

1 p

3

2

2

s ⳱ ixf/Vs, p ⳱ ixf/Vp .

p(7)

For the static case, equation (3) reduces toU(z) ⳱zD
8(1 ⳮ ␯)

r
r 2
2(1
ⳮ
2
␯
)
Ⳮ
3
dr,
f3
f

冮 冦
R

0

冢冣冧

(8)

where ␯ is the Poisson ratio. Note that since the system isaxially symmetric, the fault integrations of equations (3) and(8) are one-dimensional along the r axis in the range from0 to R. When z is directly on the fault (i.e., z ⳱ 0), theintegrals in equations (3) and (8) can be evaluated analytically by using Cauchy’s principal values, and their value ishalf of the slip (e.g., Kane, 1994),U(z ⳱ 0, x) ⳱

D
.
2

(9)

On the other hand, when z is much larger than R, equation (8) can be approximated as
U(z) ⳱ D1 ⳮ 2␯ R 2,
8(1 ⳮ ␯) z

冢冣(R K z).

(10)Note that the static displacement attenuates as the inverseof z2.
Figure 3 shows the integrands of equations (3) and (8)as a function of r using a fault model with R ⳱ 10 km andD ⳱ 1 m for the two cases of (a) z ⳱ 1 m and (b) z ⳱
10 m. The material properties of the medium are Vs ⳱ 3
km/sec, Vp ⳱ 5 km/sec, and the corresponding ␯ is 0.2188.
The real parts of the dynamic integrands of equation (3) areplotted in dotted lines for a frequency of 10 Hz, while thestatic integrands of equation (8) are shown in solid lines.
The left figure in Figure 3a shows the integrands from r ⳱
0 to 10,000 m for the case of z ⳱ 1 m. Even though wenormalized their amplitudes by their maximum values, these

Figure 3. The near-singular behaviors of the integrands of equations (3) and (8) for observation
points at (a) z ⳱ 1 m and (b) z ⳱ 10 m. The solidand dotted lines are the static and dynamic (f ⳱
10 Hz) integrands, respectively, and the dashed lineis the value for which the static integrand was subtracted from the dynamic integrand.

1158
amplitudes can barely be seen. The right drawing in Figure
3a shows the same integrand as on the left, except that the
integration range has been narrowed to the interval from 0
to 30 m. The spikelike sharp peaks of the integrand thenbecome apparent over an extremely narrow range close tor ⳱ 0 m and quickly decrease with r. This makes it clearthat in order to evaluate the integrals, close attention mustbe paid to these spikes. Otherwise, the computed amplitudeof the displacements near the fault will be in error and far
too small. By using a dense distribution of integration points,
both equations (3) and (8) give 0.500 m, which is half of theslip, as expected from equation (9). Note that the dynamic
integrands in the dotted lines are almost identical to the staticones in the solid lines. Therefore, when we subtract the static
integrand from the dynamic integrand, as shown by thedashed line, the values become nearly zero.
Similarly, Figure 3b shows the integrands for z ⳱ 10 m.
The left figure depicts the integrands from r ⳱ 0 to
10,000 m, and the right is the same integrands except for therange of 0–300 m. Again, the near singularities with sharp
peaks occur in both the static and dynamic (10 Hz) integrands in a narrow band around r ⳱ 10–30 m. The valuesof the integrations are 0.499 and 0.489 m for the static andthe dynamic cases, respectively. Similarly, since the nearsingularities are nearly the same for the static and the dynamic cases, the value of their difference becomes a smoothfunction, with the spike removed, as shown by the dashedline. Therefore, we can eliminate the near singularities of the
dynamic integrand by subtracting the static integrand.
Next, we shall evaluate the attenuation of the static displacement. Figure 4 shows the attenuation relations using
equations (8) and (10) for the same values of R and D as
before (i.e., R ⳱ 10 km and D ⳱ 1 m). When z is very close
to the fault plane, the value of displacement obtained fromequation (8) is about half that of the slip, 0.5 m, as shownin equation (9). The displacement remains nearly constantwith z up to about 1 km, around one-tenth of R. It then starts
to decrease with increasing z, on the order of 1/z. Once z
exceeds about half the radius, the displacement starts to attenuate as the order of (1/z2). Therefore, equation (10) thenbecomes a good approximation, as seen in Figure 4.
An Efficient Method for Computing Near-Fault
Ground Motions in a Layered Half-SpaceOn the basis of the results for the circular fault, we shallpropose an efficient method for carrying out the fault integration of the representation theorem for more general cases,that is, arbitrary kinematic faulting models in layered halfspaces. As shown previously, when the observation point isclose to the fault plane, the dynamic traction Green’s functions exhibit sharp peaks within the area close to the observation point. However, subtracting and adding the static
Green’s functions can eliminate the singularities of the original Green’s function. Thus, we propose the following representation theorem, instead of equation (1):

Y. Hisada and J. Bielak

Figure 4. The attenuation relation of the static
displacement using a circular fault model with R ⳱
10 km and D ⳱ 1 m in homogeneous full-space. The
displacement is nearly constant from 0 to 1 km of z
(one-tenth of the fault radius) and attenuates reversely
proportionally to the second power of z, when z exceeds about 5 km.

Uk(Y; x) ⳱冮 {T (X, Y; x) ⳮ T (X, Y)}D (X; x)dR
R

S
ik

ik

Ⳮ

i

冮 T (X, Y)D (X; x)dR,
R

S
ik

i

(11)

where Tik is the traction Green’s function of the layered halfspace at circular frequency, x, and TikS is the static traction
Green’s function of the layered half-space (x ⳱ 0). Additionand subtraction of singular integrands is common when dealing with integral evaluations or integral equations (e.g., Apsel and Luco, 1983; Colton and Kress, 1983; Hisada, 1993,
1995).
Equation (11) consists of two fault integrations. The first
integration corresponds to the dynamic Green’s function, in
which the sharp peak of the original Green’s function has
been eliminated by the subtraction of the static Green’s function. Hence, to evaluate the integral over the fault, we canlocate the integration points coarsely at regular grid nodes
on the fault plane, even when the observation point is very
close to the fault. Also, the first integration indicates that theattenuation of the slip function, Di, is governed by the dynamic Green’s functions, which consist of body and surface
waves. Thus, their amplitudes decay on the order of 1/r to1/冪r (r is the distance from the fault to the observation point).
On the other hand, the second integral in equation (11)involves the static Green’s function. To perform this fault
integration, we have to distribute the integration points
densely in the area close to the observation point in order toincorporate the contribution from the sharp peak. The CPU
time required for computing the static function, however, ismuch less than for the corresponding dynamic function. And
more importantly, since the values of the static functions
remain invariant for all frequencies, these functions need to

1159

A Theoretical Method for Computing Near-Fault Ground Motions in Layered Half-Spaces

be evaluated only once. The second integral describes theattenuation of the slip function due to the static traction of
the Green’s functions. Thus, their amplitudes would be expected to decrease rapidly, on the order of 1/r2 beyond acertain distance from the fault, as indicated in Figure 4.
While we could have chosen some alternative singular functions to remove the singularity from the dynamic Green’s
function, we have used the static Green’s functions both formathematical convenience and for the physical insight into
the fling step and rupture directivity phenomena one gains
from their use.
Efficient Methods for Calculating the Dynamic
and Static Green’s Functions When the Depth ofSource Point is Equal or Close to the Depth ofObservation Point
Next, we address another problem that arises when one
computes the Green’s functions of layered half-spaces, considering surface faulting. As mentioned earlier, when thedepth of the source is the same as or close to the depth ofthe observation point, the integrand of the wavenumber integration of equation (2) diverges with increasing wavenumbers. There are several techniques to remedy this problem,
both for the dynamic and static cases (see, e.g., Hisada,1993). We will adopt the most appropriate methods for computing the two integrals in equation (11).
One of the most efficient methods for dealing with thedivergent integrand of the dynamic Green’s function is tointroduce an asymptotic solution of the integrand at large
wavenumbers, which is equivalent to the static function.
This technique was originally introduced by Apsel and Luco
(1983), who used the integrand of the static Green’s functionof the homogeneous full-space to approximate the dynamic
function of layered half-spaces. Later, Hisada (1995) improved the method by introducing an approximate analyticalstatic solution of the layered half-space, which results in amore rapid convergence.
To evaluate the integral of the first integrand of equation(11), we can use the following wavenumber integrations instead of equation (2), since the dynamic Green’s functions
are already subtracted by the static Green’s functions:G(z,r) ⳱⬁冮 g(z)J(kr)dk⳱ 冮 g(z)J(kr)dk Ⳮ 冮 g(z)J(kr)dk,
0

kA

⬁

0

kA

⬁

⬁

kA

kA

冮 g(z)J(kr)dk ⳱ 12 冮 g(z) {H (kr)
(1)

(12)

冮 {t (X, Y; x, k) ⳮ t (X, Y; k)}dk,
0

ik

S
ik

where tik(X,Y;x) and tSik(X,Y;x) are the integrands of the dynamic and static traction Green’s functions, respectively.
Since the dynamic integrand converges exactly to the static
integrand with increasing wavenumber, the integrand of
equation (12) attenuates rapidly, even if the depth of thesource point is equal to the depth of the observation point.
Note that the expression on the right side of equation (12)is of the same form as that used by Apsel and Luco (1983)
or Hisada (1993, 1995). However, here we adopt the static(13)

where G(z,r) is an arbitrary Green’s function, z is the depthof the source or observation points, r is the horizontal distance between the source and observation points, g(z) denotes the motion-stress vectors of z (e.g., Aki and Richards,
1980; Hisada, 1993), and J(kr) is the Bessel function of thefirst kind. The corner wavenumber, kA, can be arbitrary, but
it must be larger than the wavenumbers corresponding to thebranch points and the poles in the dynamic case. For thestatic case, since there are no branch points or poles, anarbitrary small value is possible, such as kAr ⳱ 1 or 2.
The second integrand of equation (13) diverges withincreasing k when the depths of the source and observation
points are close. Thus, we shall introduce the contour deformations (Greenfield, 1995). We first replace the Bessel
function by the Hankel functions, as follows:Ⳮ H(2) (kr)}dk,

Tik(X, Y; x) ⳮ TikS (X, Y)⬁⳱Green’s functions of the actual layered system, which guarantee the fastest convergence.
On the other hand, to evaluate the static Green’s functions of a layered half-space in the second integration ofequation (11), we employ the contour deformation method,which is a well-known technique for evaluating the wavenumber integration efficiently (e.g., Levy and Mal, 1976).
Greenfield (1995) applied this technique to evaluate Green’s
function for the case in which the depth of a source point isclose or equal to that of an observation point. This method
is more appropriate than the asymptotic method for the staticcase because the approximate asymptotic solutions, such as
the static Green’s function of the homogeneous space or Hisada’s solution (Hisada, 1995), do not converge efficiently tothe integrand of the static Green’s function, due to the infinite wavelength. In order to introduce the contour deformation, we first divide the integration range into two parts,
as shown in Figure 5,

Figure 5. The contour deformation for evaluating
the second integration of equation (13), which corresponds to the static Green’s function.

(14)

1160

Y. Hisada and J. Bielak

where H(1) and H(2) are the Hankel functions of the first andsecond kinds, respectively. Next, we expand the real wavenumber k to the complex wavenumber and apply Cauchy’s
theorem. For the integration of H(1), we carry out the contour
integration along the contour ABD in Figure 5. The integral
on BD vanishes by Jordan’s lemma, as k goes to infinity.
Moreover, since H(1) decreases exponentially on the path
AB, we can terminate the integration on AB at a suitable
point B⬘ on that path. Therefore, the integration can be approximated as⬁冮 g(z)H (kr)dk ⬇ 冮
(1)

kA

AB⬘

g(z)H(1) (kr)dki . (15a)Similarly, we carry out the contour deformation for the integration of H(2) along the contour ACD in Figure 5. We
then terminate the integration range at C⬘:⬁冮 g(z)H (kr)dk ⬇ 冮
(2)

kA

AC⬘

g(z)H(2) (kr)dki . (15b)The numerical integrations in equations (15a) and (15b) canbe readily carried out, since both integrands converge
quickly to zero with increasing ki. Therefore, the sum ofequations (15a) and (15b) gives accurate values for the second integral of equation (13).
Miscellaneous Techniques for Calculating Strong
Motions Including Permanent Offset
We will now provide some additional details on thecomputation of Green’s functions, the numerical integrations, and the calculation of the near-fault ground motions,
including static offsets.
Among various available methods for computing the
dynamic and static Green’s functions of layered half-spaces,
we use the method by Hisada (1993), which is modified from
Luco and Apsel (1983). This is an appropriate method forthis study, as it is completely free from the numerical instabilities at lower and higher frequencies through the introduction of the R/T (reflection/transmission) matrix. In addition, it gives the static solution explicitly (Luco and Apsel,1983).
For the wavenumber integrations of the dynamic
Green’s functions in equation (12), which include the poles
and the branch points, we use Simpson’s and Filon’s quadratures for small and large wavenumbers, respectively(Olver, 1964). On the other hand, for the integrations of thestatic Green’s functions in equations (13) and (15), which
involve smooth integrands, we use an adaptive Newton–
Cotes quadrature.
With regard to the fault integration in equation (11), we
divide the fault plane into rectangular subfaults and use
Gaussian quadrature in each subfault. For computing the first(dynamic) integration in equation (11), we use a variable
number of Gaussian points, up to 6 ⳯ 6, to keep at least six

points per wavelength. On the other hand, for computing thesecond (static) integration in equation (11) when an observation point is close to a subfault, we divide the subfault into
smaller subfaults and redistribute Gaussian points to evaluate the near-singularity of the static Green’s function; we
repeat this process until the value of the integration converges within a prescribed tolerance.
The static offset of displacement is calculated as follows. Since the Fourier amplitude of a steplike displacement
is infinite at 0 Hz, we first compute velocities in the frequency domain, transforming them into the time domain.
Then, we calculate the corresponding displacements by integrating the velocities numerically in the time domain. To
check the accuracy of the calculated static offset, we recall
that the Fourier amplitude of velocity at 0 Hz is the same as
the value of the static offset in the time domain, since thestatic offset is the displacement as the time goes to infinity.

Numerical Examples
Strike-Slip Model with Surface Faulting in a
Homogeneous Half-SpaceWe shall first check the basic characteristics of the nearfault ground motion using a simple strike-slip model withsurface faulting (see Fig. 6). We compute the near-fault
ground motions using the proposed representation theorem(equation 11) from 0 Hz up to 10 Hz. As shown in Figure
6a, we use the 12 observation points on the free surface
along a line perpendicular to the fault plane. Point 1 is
0.1 km away from the center of the surface fault, and point
2 is 0.5 km away. Similarly, points 3–12 are 1.5–10.5 km
away from the fault at 1-km intervals. Because of the symmetry of the model and the pure strike slip, the fault-normalcomponents are zero, and the vertical components are negligible compared to the fault-parallel components. As shownin Figure 6b, the maximum slip of the fault is 1 m, including
the shallowest subfaults which break the free surface, andtapers at both edges and at the bottom of the fault. The slip
velocity function is an isosceles triangle with a 1-sec duration, as seen in Figure 6c. The rupture velocity is infinite.
We use the homogeneous half-space with the physical properties of q ⳱ 2.5 g/cm3, Vp ⳱ 5 km/sec, Vs ⳱ 3 km/sec,
Qp ⳱ 200, and Qs ⳱ 100.
Figure 7 shows the velocities (Fig. 7a) and displacements (Fig. 7b) of the fault-parallel components. Both panels
include three sets of waves on the line; the left and middle
waveforms correspond to the dynamic and static terms,
which stem from the first and second integrals in equation(11), respectively. The rightmost panels represent the sumsof the two terms. When observation points are close to thefault (e.g., points 1 and 2), we see the strong fling effects,
that is, the large pulse in velocities and the large static offsets
in displacement. These consist mainly of static terms, withlittle contribution from the dynamic terms. Since the static
terms decay rapidly with distance from the fault, the flingA Theoretical Method for Computing Near-Fault Ground Motions in Layered Half-Spaces

1161

Figure 6.

(a) Strike-slip model with surface faulting and 12 observation points;(b) the slip distribution; and (c) the slip velocity function.

Figure 7. (a) Velocities and (b) displacements of the fault-parallel components at
12 observation points in Figure 6, using the first (dynamic; left), second (static; center),and total (right) integrations of equation (11).
effects disappear for observation points away from the fault(e.g., points 11 and 12). By contrast, the dynamic terms become dominant far from the fault. Note that the static velocities are always isosceles triangles, which correspond toslip velocity function with attenuation (see Fig. 6c).
To check the attenuation relations of the dynamic, static,
and total terms with distance, Figure 8 shows the relation
between r (the distance from the fault) and the maximum
amplitudes of the velocities (Fig. 8a) and displacements(Fig. 8b). In each panel, the thin black lines, the thick graylines, and the thick black lines correspond to the dynamic,

static, and total terms, respectively. In addition, the dashedline in Figure 8b represents the Fourier amplitudes of velocity at 0 Hz, which agrees with the thick gray line, that is, thevalues of the static offset. This demonstrates the validity of
the results, as explained in the previous section. For smallr, the static terms are dominant over the dynamic ones, andthe maximum amplitudes are about half of the slip function,that is, 100 cm/sec in velocity and 50 cm in displacement.
As the distance becomes large, say, 2–5 km (one-fifth to onehalf of the fault dimension), the static terms decrease as
1/r2. By contrast, since the dynamic terms are attenuated on

1162

Y. Hisada and J. Bielak

Figure 8. Attenuation relations of the maximum amplitudes in (a) velocities and(b) displacements using the dynamic, static and total waves; 100 cm/sec in velocity
and 50 cm in displacement are the half of the slip velocity and the slip. In panel (b),since the maximum displacements of the static terms (i.e., the permanent offsets) andthe Fourier amplitudes of velocities at 0 Hz show almost same values, there are nodifferences in the figure.

the order of 1/r to 1/冪r, they become dominant over the staticones at distances greater than 2–5 km. Therefore, the fling
effects, which correspond to mostly the static terms, areprominent only in the vicinity of the fault (or asperity size).
Dip-Slip Model with Surface Faulting in a
Homogeneous Half-Space
Next, we shall check the near-fault strong motions using
a dip-slip model with surface faulting. As shown in Figure
9, the dip and rake angles are 45⬚ and 90⬚, respectively. The
other source and medium parameters are the same as for theprevious model (see Fig. 6). We locate 14 observation points
on the free surface along the line perpendicular to the fault;
points 1–7 are located on the footwall side, and points 8–14
are on the hanging-wall side. The points closest to the fault
are points 7 and 8, which are at 0.1 km away from the surface
fault trace.
Figure 10a,b shows the fault-normal and up–down components, respectively, of the velocities at the 14 observation
points. Similarly, Figure 10c,d shows the corresponding displacements. As in the previous example, the waves in theleft and middle panels in each set are the dynamic and static
terms, respectively, and the right waves are the sum of thetwo. Compared with the waves on the footwall (points 1–7),
the waves on the hanging wall (points 8–14) show large
amplitudes, especially in the up–down components. Whenobservation points are close to the fault, the static terms are
dominant over the dynamic ones, and we see strong fling
effects. As the observation points get farther from the fault,

Figure 9. Dip-slip model with surface faulting inhomogeneous half-space. The dip angle is 45⬚, andthe other source parameters are the same as Figure 6.
the static terms are quickly attenuated and the dynamic ones
become dominant.
Fling versus Directivity Effects Using Surface
and Buried Fault Models
We now examine the effects of the fling step and therupture directivity on the near-fault strong motion for strikeslip models. As shown in Figure 11a,b, we consider both asurface fault model and a buried fault model. As previous
studies have indicated, the depth of the top of the fault wouldA Theoretical Method for Computing Near-Fault Ground Motions in Layered Half-Spaces

1163

Figure 10. Velocities (a,b) and displacements (c,d) in the fault-normal (a,c) andup–down (b,d) directions at 14 observation points in Figure 9, using the dip-slip model,
and the first (dynamic; left), second (static; center), and total (right) integrations ofequation (11). The horizontal line indicates the location of the surface fault.
be one of the most important parameters controlling thenear-fault ground motions (e.g., Anderson and Luco, 1983a,
1983b; Luco and Anderson, 1983). Most source parameters
are the same as in Figure 5, except that now the rupture
velocity has a finite value of 2.5 km/sec; this gives rise todirectivity effects. The location of the hypocenter is shownin the figure, halfway down the fault. The depth of the topof the buried fault is 2 km, as shown in Figure 11b. We
locate 13 observation points at 1-km intervals on the line
parallel to the fault plane at a distance of 0.1 km. Thus,
observation points 1 and 13 are in the backward and forward
directions of the rupture front, respectively.

Figure 12 shows the results of the surface fault model,
where panels (a) and (b) are the velocities of the fault-normal
and fault-parallel components, respectively, and panels (c)and (d) are the corresponding displacements. As before, theleft, middle, and right waves in each figure are the dynamic,
static, and total contributions, respectively. In the faultnormal components, we see the growth of the directivity
pulses in the forward direction of the rupture front, which
consist mostly of the dynamic terms. The contribution of thestatic terms is minor. By contrast, we see fling waves in thefault-parallel components; these consist of the static terms
and are dominant only at the observation points next to the

1164

Y. Hisada and J. Bielak

bined effects of fling and directivity. On the other hand, thefault-normal components become dominant for observation
points further away from the surface fault, especially in theforward directivity direction.
Effects of Sedimentary Layers on Near-Fault
Ground Motions

Figure 11. A strike-slip model with (a) surface
fault and (b) buried fault in homogeneous half-space.
The source parameters are the same as Figure 6, butwith Vr ⳱ 2.5 km/sec.
fault plane (i.e., points 3–11 in Fig. 12). The contributionsof the dynamic terms are negligible.
Figure 13a,b depicts the velocities and Figure 13c,d thecorresponding displacements for the buried fault model. Notice that the scales of the amplitudes are about half of thosein Figure 12. From these figures we can see that the faultnormal components, which are mostly the dynamic waves,
exhibit a clear forward directivity effect. On the other hand,the fling waves disappear in this model. As explained before,
the fling consists mainly of static Green’s functions, which
attenuate more quickly than the dynamic ones. And, more
importantly, when an observation point is located next to asurface fault, the slip dislocation directly “flings” the adjacent ground in the slip direction. However, when an observation point is located above a buried fault, the slip dislocation cannot fling the ground, because the medium iscontinuous above the fault.
In order to further illustrate the combined effects of flingand directivity on the near-fault ground motion,Figure 14
shows the vectors of the maximum velocities on the freesurface using the surface fault (Fig. 14a) and the buried fault(Fig. 14b), as shown in Figure 11a,b. Again, the amplitude
scale for the buried fault is about half of that for the surface
fault. In both models, fault-normal components grow in theforward rupture direction (the left side of the figures). In
particular, the normal components are dominant in the vicinity of the buried fault. By contrast, we also see large amplitudes in the fault-parallel components near the surface
fault in Figure 14a. Therefore, because of the comparable
amplitude of fault-parallel and fault-normal components, thedirections of the maximum velocities are inclined with respect to the fault plane in the vicinity of the surface fault. In
particular, the fault-parallel components are larger than thenormal in the area close to the epicenter. As mentioned in
the Introduction, the direction of the maximum velocity recorded at Lucerne Valley was inclined to the fault plane (seeFig. 1b). This phenomenon probably occurred by the com-

Finally, we shall add sedimentary layers to the surface
fault model to check the effects of layering on the near-fault
ground motions. Table 1 shows the material properties ofthe layered half-space, in which we have added two sedimentary layers to the homogeneous half-space. We use thesame strike-slip model with surface faulting as that shownin Figure 11a. Bouchon (1979) examined the effects of layering using a model with a single layer and a strike-slip fault
and found that the amplification was much larger in the faultnormal component than in the parallel-fault component.
Figure 15a,b shows the three components of velocities
and displacements, respectively, at the same observation
points shown in Figure 11a. Compared with the results of
the homogeneous medium shown in the total waves of Figure 12, the fault-normal components show larger amplitudes
and longer durations, especially in the forward rupture direction. This is because the fault-normal components aremostly the dynamic terms, that is, the body and surface
waves, and they are amplified in the sedimentary layers. For
the fault-parallel components, which are represented mostly
by the static terms, the waveforms and amplitudes are nearly
the same as those of the homogeneous medium.

Concluding Remarks
In this article, we first examined a simple circular fault
model (Fig. 2) with the objective of developing efficientmethods for computing near-fault strong ground motions.
From these simple examples, we have made some important
observations. When an observation point is close to the fault
plane, the integrands of the fault integrations exhibit nearsingularities, that is, sharp peaks in the narrow area close to
the observation point, as seen in Figure 3. Therefore, in orderto obtain accurate values of the integrations without further
modification, integration points need to be distributed very
densely in the area. However, subtracting the integrand of
the static Green’s functions can eliminate the near singularity of the dynamic Green’s function. We also found that thestatic displacement (the permanent offset) attenuates rapidly
on the order of 1/r2 at distances larger than a fraction of theradius R, as shown in Figure 4.
On the basis of these observations, we have proposeda new form of the representation theorem, described byequation (11), from which one can compute the near-fault
ground motions in layered half-spaces much more efficiently
than with the original equation (equation 1). The proposed
theorem entails two fault integrations; the first consists ofthe difference between the dynamic Green’s function andthe corresponding static Green’s function, and the second
consists of the static Green’s function. In the first integration,

A Theoretical Method for Computing Near-Fault Ground Motions in Layered Half-Spaces

1165

Figure 12. Velocities (a,b) and displacements (c,d) in the fault-normal (a,c) and-parallel (b,d) directions at 13 observation points in Figure 11a, using the first (dynamic;left), second (static; center), and total (right) integrations of equation (11).

the near singularities of the dynamic Green’s functions arecompletely eliminated by the subtraction of the static
Green’s functions. Therefore, we can carry out numericallythe first integration by distributing the integration points on
regular grids that consider only the target wavelength. As
for the second fault integration, the integration points need
to be distributed very densely in the area close to the observation points in order to evaluate the near singularities. However, the required computations are small, because (1) thestatic Green’s functions are much simpler than the corresponding dynamic ones and, more importantly, (2) the computation for the static Green’s functions needs to be per-formed only once; these are then used for all frequencies. It
should be noted that the proposed representation theorem ismeant to be used only for observation points that are near
the fault. Since the near singularities of the dynamic Green’s
functions disappear as the distance between the observation
point and the fault increases, it is unnecessary to introduce
the static Green’s functions.
We also proposed efficient methods for computing the
dynamic and static Green’s functions of layered half-spacesin the case of a surface fault. Since the depths of source
points are on or close to the free surface in this case, theintegrands of the wavenumber integration in equation (2) do

1166

Y. Hisada and J. Bielak

Figure 13. Velocities (a,b) and displacements (c,d) in the fault-normal (a, e) and-parallel (b,d) directions at 13 observation points in Figure 11b, using the first (dynamic;left), second (static; center), and total (right) integrations of equation (11).
not converge. In order to carry out the wavenumber integration of the dynamic Green’s function in the first fault integration in equation (11), we use equation (12); its integrandquickly attenuates with the wavenumbers. For the static
Green’s function of the second fault integration in equation(11), we use the contour deformation method (Greenfield,
1995), as shown in equations (14) and (15). The resulting
integrands then have an asymptotic exponential decay andcan be easily integrated numerically.
Next, we computed the near-fault ground motions using
various idealized models and drew the following conclusions. First, we found that the fling step stems mainly fromthe second fault integral in equation (11), which is the convolution of the fault slip and the static traction Green’s function. Their amplitudes are dominant in the slip direction in
the vicinity of surface faults and are negligible for buried
faults. This is because the static Green’s function quickly
attenuates with distance, on the order of 1/r2. In addition,
when an observation point is located above a buried fault,
the medium has to remain continuous, and thus cannot fling.
On the other hand, we found that the directivity pulse
mainly stems from the first integral in equation (11), whichis the convolution of the fault slip and the dynamic Green’s
function. Directivity effects are dominant in the fault-normal

1167

A Theoretical Method for Computing Near-Fault Ground Motions in Layered Half-Spaces

Figure 14.

Vectors of maximum velocities on the free surface for the (a) surfaceand (b) buried faults. In the case of the strike-slip fault, the effects of “fling” and“directivity” can been seen in the components of the fault parallel and the normal,
respectively.

Table 1The Material Properties of the Layered Half-Space
Density(t/m3)

2.2
2.4
2.5

Vp
(km/sec)

Vs(km/sec)

Qp

Qs

Thickness
(km)

3.0
4.0
5.0

1.0
2.0
3.0

100
150
200

50
75
100

1
1
–

direction, especially in the forward rupture direction, notonly for the surface fault, but also for the buried fault. They
attenuate much more slowly than the fling, on the order from1/r to 1/冪r. Therefore, the fault-normal components aremuch larger than the fault-parallel components, especiallyfor the buried fault, as shown in Figure 14b. In the vicinityof surface faulting, the combined effects of directivity andfling cause the directions of maximum velocities to become
inclined, as shown in Figure 14a. This probably explains the
direction of the maximum velocity recorded at Lucerne Valley, as shown in Figure 1b.
Also, we investigated the effects of softer surface layers
using the same fault model and found that the softer layers
amplify mainly the directivity pulses, but not the fling steps.
This is because the body and surface waves of the dynamic
Green’s function are excited more strongly in the softer layers than in the static one. In actual situations, however, theexistence of softer layers may not only increase the faulting

slip, but also change the slip function itself. The results ofdynamic faulting models, laboratory tests, and/or source inversions deduced from strong motion records should be incorporated into future work to model slip functions of surface faults more appropriately.
Finally, we offer the Fortran code of this method for
public use. Please contact the first author for information on
how to obtain and use the code.

Acknowledgments
This research was partly supported by a special project of U.S. Japan
Cooperative Research for Urban Earthquake Disaster Mitigation (no.
11209201), Grant-in-Aid for Scientific Research on Priority Area (CategoryB), Ministry of Education, Culture, Sports, Science, and Technology ofJapan (MEXT), and the Earthquake and Environmental Research Center ofKogakuin University funded by the Frontier Research Promotion Programof MEXT (Y.H.). This work was also supported by the National Science
Foundation KDI program, under Grant CMS-9980063 (J.B.); the cognizant
program director is Clifford Astill. We are grateful for this support. We
obtained the strong motion records used in this study from the PEER Strong
Motion Database.

References
Abrahamson, N. (2001). Incorporating effects of near fault tectonic deformation into design ground motions, a presentation sponsored by the
Earthquake Engineering Research Institute Visiting Professional Program, hosted by the University at Buffalo, 26 October 2001, http://
mceer.buffalo.edu/outreach/pr/abrahamson.asp.

1168

Y. Hisada and J. Bielak

Figure 15. (a) Velocities and (b) displacements using the surface faulting model of
Figure 11a in the layered structure, as shown in Table 1.

Aki, K. (1968). Seismic displacement near a fault, J. Geophys. Res. 64,
321–342.
Aki, K., and P. G. Richards (1980). Quantitative Seismology, Vol. 1, W. H.
Freeman, New York.
Anderson, J. G. (1976). Motions near a shallow rupturing fault: evaluationof the effects due to the free surface, Geophys. J. 16, 575–593.
Anderson, J. G., and J. E. Luco (1983a). Parametric study of near-field
ground motions for strike-slip dislocation models, Bull. Seism. Soc.
Am. 73, 23–43.
Anderson, J. G., and J. E. Luco (1983b). Parametric study of near-field
ground motions for oblique-slip and dip-slip dislocation models, Bull.
Seism. Soc. Am. 73, 44–57.
Apsel, R. J., and J. E. Luco (1983). On the Green’s functions for a layered
half-space. II, Bull. Seism. Soc. Am. 73, 931–951.
Archuleta, R. J., and G. A. Frazier (1978). Three-dimensional numericalsimulation of dynamic faulting in a half-space, Bull. Seism. Soc. Am.
68, 541–572.
Archuleta, R. J., and S. H. Hartzell (1981). Effects of fault finiteness onnear-source ground motion, Bull. Seism. Soc. Am. 71, 939–957.
Bouchon, M. (1979). Predictability of ground displacement and velocity
near an earthquake fault: an example—the Parkfield earthquake of
1966, J. Geophys. Res. 84, 6149–6156.
Chen, X. (1995). Near-field ground motion from the Landers earthquake.
Earthquake Engineering Research Institute, Report 95-02, Caltech,
http://resolver.library.caltech.edu/caltechEERL: 1995.EERL-95-02.
Colton, D., and R. Kress (1983). Integral Equation Methods in Scattering
Theory, Wiley, New York.
Greenfield, R. J. (1995). Comment on “An efficient method for computing
Green’s functions for a layered half-space with sources and receivers
at close depths” by Y. Hisada, Bull. Seism. Soc. Am. 85, 1523–1524.
Haskell, N. A. (1969). Elastic displacements in the near-field of a propagating fault, Bull. Seism. Soc. Am. 59, 865–908.
Hisada, Y. (1993). An efficient method for computing Green’s functions
for a layered half-space with sources and receivers at close depths, I,
Bull. Seism. Soc. Am. 84, 1456–1472.
Hisada, Y. (1995). An efficient method for computing Green’s functions
for a layered half-space with sources and receivers at close depths, II,
Bull. Seism. Soc. Am. 85, 1080–1093.
Kane, J. H. (1994). Boundary Element Analysis in Engineering Mechanics,
Prentice Hall, New York.
Kawasaki, I. (1975). On the dynamic process of the Parkfield earthquakeof June 28, 1966, J. Phys. Earth 23, 127–144.

Kawasaki, I., Y. Suzuki, and R. Sato (1973). Seismic waves due to a shear
fault in a semi-infinite medium, I, J. Phys. Earth 21, 251–284.
Levy, N. A., and K. Mal (1976). Calculation of ground notion in a threedimensional model of the 1966 Parkfield earthquake, Bull. Seism. Soc.
Am. 66, 405–423.
Luco, J. E., and J. G. Anderson (1983). Steady-state response of an elastic
half-space to a moving dislocation of finite width, Bull. Seism. Soc.
Am. 73, 1–22.
Luco, J. E., and R. J. Apsel (1983). On the Green’s functions for a layered
half-space. I, Bull. Seism. Soc. Am. 73, 909–929.
Oglesby, D. D., and S. M. Day (2001). Fault geometry and the dynamics
of the 1999 Chi-Chi (Taiwan) earthquake, Bull. Seism. Soc. Am. 91,
1099–1111.
Olver, F. W. J. (1964). Bessel functions of integer order, in Handbook ofMathematical Functions with Formulas, Graphs, and Mathematical
Tables, Applied Mathematics Series, Vol. 55, M. Abramowitz andI. A. Stegun (Editors), U.S. Department of Commerce, National Bureau of Standards, Washington, D.C.
Ramancharla, P. K., and K. Meguro (2001). Numerical modeling of dipslip faults for studying ground surface deformation, Bull. Earthquake
Resistant Structure Research Center 34, 107–112.
Shin, T. C., and T. L. Teng (2001). An overview of the 1999 Chi-Chi,
Taiwan, earthquake, Bull. Seism. Soc. Am. 91, 895–913.
Somerville, P. G., N. F. Smith, R. W. Graves, and N. A. Abrahamson(1997). Modification of empirical strong ground motion attenuation
relations to include the amplitude and duration effects of rupture directivity, Seism. Res. Lett. 68, 199–222.
Department of Architecture
Kogakuin University
Nishi-Shinjuku 1-24-2
Tokyo 163-8677, Japan
hisada@cc.kogakuin.ac.jp(Y.H.)
Department of Civil and Environmental Engineering
Carnegie Mellon University
Pittsburgh, Pennsylvania 15213
jbielak@cmu.edu(J.B.)
Manuscript received 1 August 2002.