UCLA
UCLA Previously Published Works
Title
Use of Exact Solutions of Wave Propagation Problems to Guide Implementation of Nonlinear
Seismic Ground Response Analysis Procedures
Permalink
https://escholarship.org/uc/item/12f5332m
Journal
J. Geotech. & Geoenv. Engrg., 133(11)
Authors
Kwok, Annie O.L.
Stewart, Jonathan P
Hashash, Youssef
et al.
Publication Date
2007
Peer reviewed
eScholarship.org Powered by the California Digital Library
University of California
Use of Exact Solutions of Wave Propagation Problems
to Guide Implementation of Nonlinear Seismic Ground
Response Analysis Procedures
Annie O. L. Kwok, M.ASCE
1
; Jonathan P. Stewart, M.ASCE
2
; Youssef M. A. Hashash, M.ASCE
3
;
Neven Matasovic, M.ASCE
4
; Robert Pyke, M.ASCE
5
; Zhiliang Wang, M.ASCE
6
; and
Zhaohui Yang, M.ASCE
7
Abstract: One-dimensional nonlinear ground response analyses provide a more accurate characterization of the true nonlinear soil
behavior than equivalent-linear procedures, but the application of nonlinear codes in practice has been limited, which results in part from
poorly documented and unclear parameter selection and code usage protocols. In this article, exact 共linear frequency-domain兲 solutions for
body wave propagation through an elastic medium are used to establish guidelines for two issues that have long been a source of
confusion for users of nonlinear codes. The first issue concerns the specification of input motion as “outcropping” 共i.e., equivalent
free-surface motions兲 versus “within” 共i.e., motions occurring at depth within a site profile兲. When the input motion is recorded at the
ground surface 共e.g., at a rock site兲, the full outcropping 共rock兲 motion should be used along with an elastic base having a stiffness
appropriate for the underlying rock. The second issue concerns the specification of viscous damping 共used in most nonlinear codes兲 or
small-strain hysteretic damping 共used by one code considered herein兲, either of which is needed for a stable solution at small strains. For
a viscous damping formulation, critical issues include the target value of the viscous damping ratio and the frequencies for which the
viscous damping produced by the model matches the target. For codes that allow the use of “full” Rayleigh damping 共which has two target
frequencies兲, the target damping ratio should be the small-strain material damping, and the target frequencies should be established
through a process by which linear time domain and frequency domain solutions are matched. As a first approximation, the first-mode site
frequency and five times that frequency can be used. For codes with different damping models, alternative recommendations are
developed.
DOI: 10.1061/共ASCE兲1090-0241共2007兲133:11共1385兲
CE Database subject headings: Earthquakes; Ground motion; Wave propagation; Seismic effects; Damping
.
Introduction
Nonlinear seismic ground response analysis is seldom used in
practice by nonexpert users because parameter selection and code
usage protocols are often unclear and poorly documented, the
effects of parametric variability on the analysis results are un-
known, and the benefits of nonlinear analyses relative to
equivalent-linear analyses are often unquantified. This article pre-
sents initial results of a broad study intended to resolve these
issues so as to encourage appropriate applications of one-
dimensional 共1D兲 nonlinear seismic ground response analysis
codes in engineering practice. The goals of the project are to
provide clear and well documented code usage protocols and to
verify the codes over a wide range of strain levels.
This paper considers five leading nonlinear seismic ground
response analysis codes: DEEPSOIL 共Hashash and Park 2001,
2002; Park and Hashash 2004; www.uiuc.edu/⬃deepsoil兲,
D-MOD_2 共Matasovic 2006兲, a ground response module in the
OpenSees simulation platform 共Ragheb 1994; Parra 1996;
Yang 2000; McKenna and Fenves 2001; opensees.berkeley.edu兲,
SUMDES 共Li et al. 1992兲 and TESS 共Pyke 2000兲. The study
focuses on two issues related to the application of nonlinear codes
that can be resolved by comparing the results of such analyses
共under linear condition兲 to known theoretical solutions. The first
issue concerns the specification of input motions as “outcropping”
共i.e., equivalent free-surface motions兲 versus “within” 共i.e., mo-
1
Project Engineer, Praad Geotechnical Inc., 5465 South Centinela
Ave., Los Angeles, CA 90066-6942.
2
Professor and Vice Chair, Dept. of Civil and
Environmental Engineering, Univ. of California, Los Angeles 5731
Boelter Hall, Los Angeles, CA 90095 共corresponding author兲. E-mail:
jstewart@seas.ucla.edu
3
Associate Professor, Dept. of Civil and Environmental Engineering,
RM 2230C NCEL, MC-250, Univ. of Illinois at Urbana-Champaign, 205
N. Mathews Ave., Urbana, IL 61801.
4
Associate, GeoSyntec Consultants, 2100 Main St., Ste. 150, Hunting-
ton Beach, CA 92648.
5
Consulting Engineer, 1076 Carol Lane, No. 136, Lafayette, CA
94549.
6
Senior Engineer, Geomatrix Consultants Inc., 2101 Webster St., 12th
Floor, Oakland, CA 94612.
7
Engineer, URS Corporation, 1333 Broadway, Suite 800, Oakland,
CA 94612.
Note. Discussion open until April 1, 2008. Separate discussions must
be submitted for individual papers. To extend the closing date by one
month, a written request must be filed with the ASCE Managing Editor.
The manuscript for this paper was submitted for review and possible
publication on July 11, 2006; approved on February 6, 2007. This paper is
part of the Journal of Geotechnical and Geoenvironmental Engineer-
ing, Vol. 133, No. 11, November 1, 2007. ©ASCE, ISSN 1090-0241/
2007/11-1385–1398/$25.00.
JOURNAL OF GEOTECHNICAL AND GEOENVIRONMENTAL ENGINEERING © ASCE / NOVEMBER 2007 / 1385
tion occurring at depth within a site profile兲. The second issue
concerns the specification of the damping that occurs within a soil
element at small strains, which is either accomplished using vis-
cous damping or unload-reload rules that produce nonzero small-
strain hysteretic damping.
This paper begins with a brief review of frequency-domain
and time-domain ground response analysis procedures. This is
followed by sections describing verification studies addressing the
issues of input motion specification and modeling of small-strain
damping.
One-Dimensional Ground Response
Analysis Procedures
In 1D seismic ground response analyses, soil deposits are as-
sumed to be horizontally layered over a uniform half-space. The
incident wave is assumed to consist of vertically propagating
shear waves. The response of a soil deposit to the incident motion
can be modeled in the frequency or time domains, as described
below.
Frequency-Domain Analysis
Frequency domain analyses are based on a closed form solution
of the wave equation for shear wave propagation through a lay-
ered continuous medium, with each layer i having a specified
density
i
, shear modulus G
i
, and hysteretic damping 
i
. The
solution was presented by Roesset and Whitman 共1969兲, Lysmer
et al. 共1971兲, Schnabel et al. 共1972兲, and is also described in detail
by Kramer 共1996兲. In these frequency-domain methods, a control
motion of frequency is specified at any layer j in the system.
An exact solution of the system response can be expressed as a
transfer function relating the sinusoidal displacement amplitude in
any arbitrary layer i to the amplitude in layer j
F
ij
=
a
i
共兲 + b
i
共兲
a
j
共兲 + b
j
共兲
共1兲
where F
ij
⫽amplitude of transfer function between layers i and j;
a
i
and a
j
⫽normalized amplitudes of upward propagating waves
in layers i and j; and b
i
and b
j
⫽normalized amplitudes of
downward propagating waves in layers i and j. The normaliza-
tion of the wave amplitudes is generally taken relative to the
amplitude in layer 1, for which a
1
=b
1
due to perfect wave reflec-
tion at the free surface. The normalized amplitudes a
i
, a
j
, b
i
,
and b
j
can be computed from a closed-form solution of the
wave equation, and depend only on profile characteristics 共i.e.,
material properties , G, and  for each layer and individual layer
thicknesses兲.
The frequency domain solution operates by modifying, rela-
tive to the control motion, the wave amplitudes in any layer i for
which results are required. These analyses are repeated across all
the discrete frequencies for which a broadband control motion is
sampled, using the fast Fourier transform. Once amplitudes a
i
and
b
i
have been computed for a given layer at all those frequencies,
time-domain displacement histories of layer i can be calculated
by an inverse Fourier transformation.
Control motions for use in frequency domain analyses are
most often recorded at the ground surface, and are referred to as
“outcropping.” As perfect wave reflection occurs at the ground
surface, incident and reflected wave amplitudes are identical,
and hence outcropping motions have double the amplitude of
incident waves alone. Consider the example in Fig. 1. Rock layer
n occurs at the base of a soil column in Case 1 and as outcropping
rock in Case 2. In the outcropping rock case, incident and re-
flected waves are equivalent 共a
n
*
=b
n
*
兲. The incident waves are
identical in both cases 共a
n
*
=a
n
兲, assuming equal rock moduli, but
the reflected waves differ 共b
n
*
⫽b
n
兲 because some of the incident
wave transmits into the soil 共nonperfect reflection兲 for Case 1,
whereas perfect reflection occurs in Case 2. The motion at the
base of the soil column in Case 1 共referred to as a “within” mo-
tion兲 can be evaluated from the outcropping motion using the
transfer function
F
nn
* =
u
n
u
n
*
=
a
n
共兲 + b
n
共兲
2a
n
共兲
共2兲
As with any other transfer function, F
nn
* can be readily computed
for any frequency and depends only on profile characteristics.
Accordingly, through the use of Eq. 共2兲, the within motion can be
calculated for a given outcropping motion. The base-of-profile
共within兲 motion can in turn be used to calculate motions at any
other layer per Eq. 共1兲.
The application of Eq. 共2兲 results in a within motion that is
reduced from an outcropping motion at the site 共modal兲 frequen-
cies. Consider, for example, a single soil layer with thickness⫽
30 m, V
s
=300 m/s 关giving a fundamental mode site frequency of
共f
s
=300 m/s兲/共4 ⫻ 30 m兲=2.5 Hz兴 overlying a half-space with
shear wave velocity V
s−H
. The results of the within/outcropping
calculation 共i.e., Eq. 共2兲兲 are shown in Fig. 2共a兲 for various values
of equivalent viscous damping ratio 共equal damping values are
applied in both the soil layer and half-space兲 with V
s−H
=2V
s
and
in Fig. 2共b兲 for zero damping and various levels of velocity con-
trast 共V
s
/V
s−H
兲. As shown in Fig. 2共a and b兲, the transfer function
amplitude 共within/outcropping兲 drops below unity near the site
frequencies, with the amplitudes at site frequencies decreasing
only with decreasing amounts of equivalent viscous damping. At
frequencies between the site frequencies, amplitudes decrease
both with increasing damping and with decreasing velocity
contrast.
At zero damping the transfer function amplitude goes to zero
at site frequencies. To understand this phenomenon, consider that
共1兲 control motion and response are in phase in this case because
of the lack of damping, and 共2兲 the site frequencies correspond to
2n + 1 quarter-wave lengths, where n = 0, 1, 2, etc. 共zero and posi-
tive integers兲. As shown in Fig. 2共c兲, at a depth below the surface
of 2n + 1 quarter-wave lengths, the wave amplitude is zero 共i.e.,
there is a “node” in the response at that depth兲, which in turn must
produce a zero transfer function amplitude Fig. 2共c兲 shows mode
shapes for the first and third modes, i.e., n =0 and 1. Additionally,
as shown in Fig. 2共c兲共lower right frame兲, as damping increases,
the input and response are increasingly out of phase, and there are
no true nodes in the site response.
Fig. 1. Incident and reflected waves in base rock layer for case
of soil overlying rock and outcropping rock 共amplitudes shown are
relative to unit amplitude in Case 1 surface layer兲
1386 / JOURNAL OF GEOTECHNICAL AND GEOENVIRONMENTAL ENGINEERING © ASCE / NOVEMBER 2007
The trends shown in Figs. 2共a and b兲 at frequencies between
the site frequencies can be explained as follows: 共1兲 The decrease
of within motion amplitude with increasing damping results from
a reduction of reflected energy from the ground surface as damp-
ing increases, thus reducing the amplitude of within motions 共that
are the sum of incident and reflected waves兲; and 共2兲 the decrease
of within motion amplitude with decreasing V
s−H
results from
increased transmission of reflected 共downward propagating兲
waves from the surface into the halfspace 共i.e., less reflection兲,
which causes energy loss from the system.
Time-Domain Analysis
The principal limitation of traditional frequency domain analysis
methods is the assumption of constant soil properties 共G and 兲
over the duration of earthquake shaking. Time-domain analysis
methods allow soil properties within a given layer to change with
time as the strains in that layer change. Modified frequency-
domain methods have also been developed 共Kausel and Assimaki
2002; Assimaki and Kausel 2002兲 in which soil properties in
individual layers are adjusted on a frequency-to-frequency basis
to account for the strong variation of shear strain amplitude with
frequency. Since the frequencies present in a ground motion
record vary with time, this can provide a reasonable approxima-
tion of the results that would be obtained from a truly nonlinear,
time-stepping procedure. Nonetheless, the present focus is on
true, time-stepping procedures.
The method of analysis employed in time-stepping procedures
can, in some respects, be compared to the analysis of a structural
response to input ground motion 共Clough and Penzien 1993;
Chopra 2000兲. Like a structure, the layered soil column is ideal-
ized either as a multiple degree of freedom lumped mass system
关Fig. 3共a兲兴 or a continuum discretized into finite elements with
distributed mass 关Fig. 3共b兲兴. Whereas frequency-domain methods
are derived from the solution of the wave equation with specified
boundary conditions, time-domain methods solve a system of
coupled equations that are assembled from the equation of mo-
tion. The system is represented by a series of lumped masses or
discretized into elements with appropriate boundary conditions.
Table 1 summarizes the manner in which mass is distributed and
nonlinear behavior is simulated for the five nonlinear codes con-
sidered here.
The system of coupled equations is discretized temporally and
a time-stepping scheme such as the Newmark  method is em-
ployed to solve the system of equations and to obtain the response
at each time step. TESS utilizes an explicit finite-difference solu-
tion of the wave propagation problem that is the same as the
solution scheme used in FLAC developed by HCItasca. Unlike in
Fig. 2. Ratio of within to outcropping amplitudes for: 共a兲 various equivalent viscous damping ratios; 共b兲 various base layer velocities 共V
s−H
兲; and
共c兲 mode shapes for various conditions
JOURNAL OF GEOTECHNICAL AND GEOENVIRONMENTAL ENGINEERING © ASCE / NOVEMBER 2007 / 1387
frequency-domain analysis where the control motion could be
specified anywhere within the soil column, in time-domain analy-
sis, the control motion must be specified at the bottom of the
system of lumped masses or finite elements.
Specification of Input Motion
There has been confusion regarding the nature of the input motion
that should be specified for time-domain analyses at the base of
the profile. Consider the common case where the motion that is to
be applied was recorded at the surface of a rock site 共outcrop
motion兲. One school of thought that has been applied in practice
for many years is that the outcropping motion should be
converted to a within motion using frequency-domain analysis
关e.g., Eq. 共2兲兴, and that this within motion should then be speci-
fied for use at the base of the site profile for application in time-
domain analysis. Most users of this approach were aware that
the layer properties used in the outcropping-to-within conversion
were a potentially crude approximation to the actual nonlinear
soil properties. The approximation was accepted, however, due
to the lack of a practical alternative for obtaining within motions.
The second school of thought is that the outcropping rock motion
should be applied directly at the base of the site profile without
modification. Normally, this direct use of the outcropping motion
is accompanied by the use of a compliant base in the site profile
共the base stiffness being compatible with the character of the
underlying rock兲, which allows some of the energy in the vi-
brating soil deposit to radiate down into the halfspace 共Joyner
and Chen 1975兲. Rigid base options are also available in all
time-domain codes, but are seldom used because the condit-
ions under which the rigid base should be applied are poorly
understood.
To evaluate which of the two above approaches is correct,
time-domain analyses with elastic material properties are exer-
cised, for which frequency-domain analyses provide an exact so-
lution. This can be investigated using linear analyses because the
underlying issue involves the differences in linear wave propaga-
tion modeling with frequency-domain and time-domain analyses.
Consider, for example, a single soil layer with thickness⫽30 m,
shear wave velocity V
s
=300 m/s 共site frequency⫽2.5 Hz兲 that
overlies an elastic half-space with V
s−H
=2V
s
=600 m/s. Equiva-
lent viscous damping is assumed constant at 5%. A control motion
is selected to represent an extreme scenario with respect to the
variability between outcropping and within, which is a sine wave
at the site frequency. As shown in Fig. 4共c兲, the particular motion
selected has a frequency of 2.5 Hz, 12 cycles of shaking, and
cosine tapers at the beginning and end of the signal with a four-
cycle taper duration 共the tapers have the shape of half a cosine
wavelength兲. The control motion is specified for an outcropping
condition. A large suppression of the within motion relative to the
outcropping motion would be expected for this signal 共e.g., as
suggested by Fig. 2兲.
A frequency domain solution is exact because the material
properties are elastic 共i.e., strain invariant兲. The frequency domain
calculations are performed with the computer program SHAKE04
共Youngs 2004兲, which is a modified version of the original
SHAKE program 共Schnabel et al. 1972兲. Both the within motion
and the motion at the surface of the soil layer are calculated, with
the results shown in Figs. 4共a and b兲 with the solid black lines.
Linear time-domain analyses are performed for this site using
the “nonlinear” codes listed in Table 1 共the codes are imple-
mented with linear backbone curves兲. Four combinations of con-
trol motion and base condition are considered:
1. Outcropping motion 关Fig. 4共c兲兴 with elastic base 共V
s−H
=600 m/s兲.
2. Within motion 关which is extracted from frequency-domain
analysis; see Fig. 4共b兲兴 with elastic base.
3. Outcropping motion with rigid base 共V
s−H
=30,000 m/s or
select the “rigid base” option in nonlinear code, if available兲.
4. Within motion with rigid base.
The results in Fig. 4共a兲 show that the surface acceleration histo-
ries for Cases 共1兲 and 共4兲 match the known solution from
frequency-domain analysis. Using the within motion with an elas-
tic base 共Case 2兲 underestimates the surface motions, while using
the outcropping motion with a rigid base 共Case 3兲 overestimates
the surface motions.
Based on the above, our recommendations are as follows: 共i兲
For the common case in which the control motion is recorded as
outcropping, the motion should be applied without modification
for time-domain analyses with an elastic base; and 共ii兲 if time-
domain analyses are to be used to simulate the response of a
Table 1. Mass Representation and Constitutive Models Used in
Nonlinear Codes
Nonlinear
code
Mass
representation Constitutive model
D-MOD_2 Lumped mass MKZ 共Matasovic and Vucetic 1993兲
DEEPSOIL Lumped mass Extended MKZ 共Hashash and Park 2001兲
OpenSees Distributed
mass
Multiyield surface plasticity
共Ragheb 1994; Parra 1996; Yang 2000兲
SUMDES Distributed
mass
Bounding surface plasticity
共Wang 1990兲 and other models
TESS Distributed
mass
HDCP 共EPRI 1993兲
Fig. 3. 共a兲 Lumped mass system; 共b兲 distributed mass system
1388 / JOURNAL OF GEOTECHNICAL AND GEOENVIRONMENTAL ENGINEERING © ASCE / NOVEMBER 2007
