Bulletin of the Seismological Society of America, 91, 2, pp. 381–396, April 2001

Nonlinear Seismic Effects in Soils: Numerical Simulation and Studyby Olga V. Pavlenko

Abstract To study the effects of nonlinearity in the seismic response of soils, anumerical simulation of the propagation of vertically incident seismic waves in horizontal soil layers were performed. Shear noiselike and monochromatic seismic
waves of various intensities were used as input signals. The behavior of soils was
described by a nonlinear hysteretic model. To extract and study nonlinear components in the ground response, the nonlinear system identification method and analysis
of higher-order spectra of oscillations on the surface were applied. Even for weak
input signals, the response of the simulated soils contained a noticeable nonlinear
component. An increase in the intensity of input signals led to increasing distortionsof propagating signals, due to the generation and growth of combination-frequency
harmonics. The results show that odd types of nonlinearity are most typical for soils,
such as cubic and fifth-order nonlinearities, causing generation of the third and fifth
higher harmonics of main frequencies of input signals. Nonlinearities of even types,
such as quadratic, fourth-order, and sixth-order, concerned with asymmetry, or skewness, of oscillations (i.e., quasi-static deformations of the surface) are usually weak,
except some special cases, in which a stress-strain relationship of a soil can be represented by functions with noticeable even components. A weak nonlinearity results
in an increase in high-frequency components, due to the generation of higher harmonics. In cases of strong nonlinearity, in which a decrease in amplification and inshear moduli become noticeable, changes in spectra of propagating signals achieve
their maximum. As a result, input signals with arbitrary spectra are transformed into
output signals with spectra of the type of E(f )  f k, where k depends on the propertiesof the medium.
Introduction
The nonlinear behavior of sedimentary soils in strong
ground motions remains an unsatisfactorily explored field.
In strong earthquakes, we can observe sharp changes in spectral composition and amplitudes of oscillations dependingon ground conditions, soil consolidation and sagging, soil
loosening and lifting, liquefaction of water-saturated soils,
residual deformations, and so on. Overall, the effect is rather
complicated and depends on the composition and thicknessof soft deposits, on their saturation with water, and on thelevel of underground water, as well as on the magnitude andfrequency content of an earthquake.
At the beginning of the 1970s, laboratory experiments
had been performed to study the behavior of soft sediments
under loading (Hardin and Drnevich, 1972), and it was foundthat the stress-strain relationship was hysteretic (Fig. 1a) andthat the behavior of sediments was elastic plastic or elastic
viscous plastic. Similar relationships were obtained fromfield experiments with small explosions, in which stress andstrain measurements in loam and loess above the level of theunderground water were performed simultaneously (Fig. 1b)(Vasil’ev et al., 1969, 1977). Downhole vertical arrays pro-vided further field evidence of nonlinear dynamic soil behavior during earthquakes. Ground motions recorded at the
Lotung large-scale seismic test site inferred reduction inshear modulus with increasing effective shear strain, which
is found to compare satisfactorily with those obtained
through laboratory tests. Nonlinear stress-strain dependencies were obtained directly from the free-field downhole accelerations (Elgamal et al., 1995; Zeghal et al., 1995; Changet al., 1996).
For a long time, there was a substantial difference in the
viewpoints of seismologists and geotechnical engineers on
the role of nonlinear wave effects in seismic fields: engineers
have been convinced of strong nonlinearities in the behaviorof soft soils (at strains smaller than 10ⳮ4), while seismologists interpreted strong motion data up to a 0.3–0.4g levelwith a linear viscoelastic behavior. This discrepancy progressively disappeared with the accumulation of representative experimental data on strong ground motion. Noticeable manifestations of soil nonlinearity were observed inrecent earthquakes (Caillot and Bard, 1990; Aki and Irikura,
1991; Darragh and Shakal, 1991; Beresnev et al., 1995;
381

382

O. V. Pavlenko

Figure 1.

(a) An initial loading curve of simple
shear and hysteresis loops for the elastic-plastic
model of soil behavior. Stress and strain are scaled sothat the maximum stress and the low-strain modulus
are unity, that is, stress is normalized by multiplying
by 1/smax, where smax is the shear stress in failure, andstrain is normalized by multiplying by Gmax /smax,
where Gmax is the low-strain modulus (derived from
Joyner and Chen, 1975). (b) Typical compression and(c) shear stress-strain curves for dry loam in (1) loading and (2) unloading (the natural humidity of 10–
15%, Vp ⳱ 300 ⳰ 600 m/sec, Vs ⳱ 200 ⳰ 250 m/
sec). Field experiments were performed with small
explosions, and the relationships were calculated
from measurements of the main stress, rii(t), andstrain, eii(t), components in soils near the surface.
Normal stresses were measured by small rigid membrane transmitters. Strains were measured by borehole strainmeters with bases of 250–300 mm. Areasof the circles are proportional to the numbers of measurements (Vasil’ev et al., 1977).

Aguirre and Irikura, 1997), and statistical analysis confirmed
the tendency toward increasing nonlinearity with increasing
intensity of oscillations (Kamiyama, 1992). To describe thenonlinear behavior of soils, various models were constructedwith elements of elasticity, plasticity, and viscosity. The conclusion was made that the direction to be followed to bridge
the gap between seismologists and geotechnical engineers isthe detailed analysis allowed by strong motion data fromdense seismic arrays (Bard and Pitilakis, 1995).
In Russia, the nonlinear wave effects in seismic fields
have been studied since the 1960s and 1970s (Nikolaev,
1967; Vasil’ev et al., 1969, 1977). First, the purpose was toevaluate nonlinear corrections to the physical characteristics(energy, spectral composition, etc.) of powerful sources ofseismic waves (seismic sources, explosions, seismic vibrators), and later to develop new types of receiving and transmitting antenna systems (Gushchin and Shalashov, 1981).
The following possible practical applications were proposed
based on the determination of nonlinear parameters of sedimentary and crystalline fractured rocks: new methods ofseismic microzonation, new sensitive methods of seismic
monitoring (based on theoretical and experimental facts that
nonlinear parameters are more sensitive to changes in the
stress state of the medium than linear parameters), and methods for seismic prospecting (Nikolaev, 1988; Groshkov etal., 1990). Nonlinear seismology developed as a continuation of nonlinear acoustics to seismic frequency bands, andmultiple experiments with a vibrator seismic source were
performed, aimed at detecting, in seismic fields, the nonlinear effects that had been already found in acoustical fields.
Interaction of seismic waves, generation and growth of combination-frequency harmonics, and other nonlinear effects,
as well as the influence of dispersion, absorption, and nonlinearity on seismic wave fields, were studied (e.g., Aleshin
et al., 1981; Gushchin and Shalashov, 1981; Nikolaev et al.,1995).
The nonlinear phenomena observed in sedimentary soils
in strong motions are due to the fact that Hooke’s law does
not hold for sediments, the stress-strain relationship is nonlinear, that is, its loading and unloading parts can be represented by the power seriesr(e) ⳱ a1e Ⳮ a2e2 Ⳮ a3e3 Ⳮ a4e4 Ⳮ a5e5 Ⳮ . . . ,(1)where r is stress, e is strain, and a1, a2, a3, a4, a5 are constantcoefficients. The first term in the series a1e describes thelinear dependence (Hooke’s law); the others are nonlinear
corrections. The nonlinearity of the stress-strain relationshipleads to distortions in the shapes of propagating seismic
waves; the waves acquire new spectral components. In otherwords, the waves interact with themselves and with eachother in a nonlinear medium, producing new components intheir spectra. The quadratic term a2e2 defines the interactionof two waves (for example, at frequencies x1 and x2) generating combination-frequency harmonics at frequencies
2x1, 2x2, x1 Ⳳ x2 and the constant component at zero fre-

383

Nonlinear Seismic Effects in Soils: Numerical Simulation and Study

quency. The cubic term a3e3 defines the interaction of threewaves at frequencies x1, x2, and x3 generating combinationfrequency harmonics at 3x1, 3x2, 3x3, x1 Ⳮ x2 Ⳮ x3,
x1 Ⳮ x2 ⳮ x3, x1 Ⳮ 2x2, 2x1 ⳮ x2, and so forth.
To account for the medium nonlinearity in problems ofwave propagation, nonlinear corrections are introduced into
wave equations, according to formula (1). By analogy withnonlinear acoustics, it was assumed that the most significant
correction to Hooke’s law is the quadratic term a2e2, andfield experiments performed in Russia were aimed at detecting the quadratic nonlinearity of soils (e.g., Aleshin etal., 1981; Gushchin and Shalashov, 1981; Groshkov et al.,
1990; Nikolaev et al., 1995). The assumption of the predominance of quadratic nonlinearity is usually made in studying
nonlinear seismic effects (Lund, 1983; Parker, 1988; Dimitrui, 1990). This article shows that soils possess mostly odd
types of nonlinearity, that is, the cubic correction to Hooke’s
law is most sufficient, whereas the quadratic and other evenorder terms are less important; however, in some cases evenorder terms become comparable with odd-order terms. It
should be mentioned that, though field experiments were
aimed at detecting even-order higher harmonics, the results
often showed a more noticeable increase in harmonics of
odd orders (Dimitriu, 1990; Nikolaev et al., 1995).
Though even and odd types of nonlinearity have muchin common, they lead to different nonlinear effects. Frequencies of higher harmonics, shapes of seismic solitary
waves, and other nonlinear phenomena depend on the typeof soil nonlinearity; only even types of nonlinearity are connected with the constant component of the seismic wave
field. Preliminary study is required to determine what corrections to the wave equations are most important in one
case or another.
In this work, to determine the types and quantitative
characteristics of soil nonlinearity, the method of nonlinear
system identification was applied, which was successfully
used in nonlinear system analysis in many fields of science.
A soil with a typical stress-strain relationship was studied to
discover general regularities in nonlinear seismic phenomena. This work will be continued by investigation of thenonlinearity of soils in various geotechnical conditions, including their liquefied states. Such work using vertical arrayrecords of the 1995 Kobe earthquake is now in preparation.
The issue concerned with the types of soil nonlinearity
is closely connected to the problem of generation of higher
harmonics of propagating seismic waves. Generation ofhigher harmonics has not yet been fully explored. Field experiments with a vibrator source were performed without
preliminary modeling, accounting for real soil properties,
whereas nonlinear models of soil behavior were only applied
to calculations of earthquakes and were not used for studying
basic nonlinear seismic phenomena. Therefore, another goalof this study is to investigate higher harmonic generation insoils by means of numerical simulation. It is shown that
nonlinearity of soil behavior can lead to a certain shape ofspectra of propagating signals, in which spectra take theform of E(f ) ⳱ f ⳮk. This is in agreement with the common
theory of wave interactions in a nonlinear medium.

Nonlinear System Identification and Higher-Order
Spectra and Their Application to the Study ofNonlinearity of Subsurface Soils
The nonlinearity of the seismic response of subsurface
soils is determined by their nonlinear stress-strain characteristics, which play the parts of the transfer functions anddetermine the transformation of input seismic signals into
the ground response. Figure 1a, b, and c shows that only the
beginning part of this curve, the domain of small deformations, is approximately linear. Therefore, only weak seismic
waves keep their shapes in propagating through such media,
whereas intense waves suffer distortions. To estimate aground seismic response at a given point, we should determine the linear and nonlinear domains in the response andcompare them with a whole set of possible seismic actions.
We should obtain a quantitative description of the groundresponse to any seismic action. In system analysis this procedure is called nonlinear system identification. The problem
of nonlinear system identification is finding such a mathematical description (model) of a system that its seismic response coincides with the response of the real physical system (Marmarelis and Marmarelis, 1978).
To obtain such a description, let us represent our system
as a black box with an input and an output. We look for a
relation between the input and output signals. In a common
case, under the conditions of stationarity, analyticity, and afinite memory of a system (these conditions are usually satisfied for real physical systems), an output signal can be
represented as the Volterra series, that is, a sum of multiple
integrals of an input signal (Marmarelis and Marmarelis,1978):
y(t) ⳱ k0 Ⳮ冮 k (s)x(tⳮs)ds Ⳮ
1

0

 

冮 冮 k (s ,s )x(tⳮs )x(tⳮs )ds ds Ⳮ
2

0

1

2

1

2

1

2

(2)

0

  

Ⳮ

冮 冮 冮 k (s ,s ,s )x(tⳮs )x(tⳮs )x(tⳮs )ds ds ds Ⳮ . . . ,
3

0

0

1

2

3

1

2

3

1

2

3

0

where x(t) is an input signal, y(t) is the output signal, t is
time, s, s1, s2, s3 are time delays, and k0, k1(s), k2(s1, s2),
k3(s1, s2, s3) are the zero-order, first-order, second-order, andthird-order Volterra kernels of the system, which are symmetric functions.
If a system is linear,k2(s1, s2) ⳱ k3(s1, s2, s3) ⳱ . . . ⳱ 0,

384

O. V. Pavlenko

and knowledge of the first term of the Volterra series (i.e.,knowledge of the first-order kernel k1(s), which is the impulse characteristic of the system) is sufficient to describe
the system. To describe a nonlinear system, several terms
are necessary. The majority of physical systems can be satisfactorily described by the first few terms because the series
converges rather quickly. This is stipulated both by thesmoothness of the functionals and by a linearizing effect ofnoise (noise accompanying measurements as well as noiseoccurring inside the system).
The Wiener series is widely used along with the Volterra series. It consists of functionals, which are linear combinations of functionals of the Volterra type. Terms in the
Wiener series are orthogonal with respect to the input signal
in the form of the Gaussian white noise, so we can determine
the Wiener kernels independently from each other. From theother side, the Gaussian white noise as an input signal allows
the most effective testing of the system because the Gaussian
white noise is a very rich signal: it contains components ofall frequencies and amplitudes. Thus, the determination of
the Wiener kernels of a system is a good method for its
identification (Marmarelis and Marmarelis, 1978). In thiscase, the dependence of the system response upon an input
signal can be written in the form:y(t) ⳱

兺 Gm[hm(s1, . . . , sm); x(t⬘), t⬘ ⱕ t],

(3)

m⳱0

where Gm are orthogonal functionals, if x(t) is the Gaussian
white noise with a zero mean, {hm(s1, . . . , sm)} is the sequence of the Wiener kernels, and s1, . . . , sm are time delays.
Each kernel is a symmetric function. The first four Wiener
functionals are represented in the form:
G0[h0; x(t)] ⳱ h0 ,(4)

h0 ⳱ E[y(t)].

(8)

The first-order kernel is expressed through the cross-correlation function of the input white noise and the response of
the system:
h1(s) ⳱ (1/P)E[y(t)x(t ⳮ s)].

(9)

The second-order kernel is defined ash2(s1, s2) ⳱ (1/2P2)E[y1(t)x(t ⳮ s1)x(t ⳮ s2)],

(10)

where



G1[h1; x(t)] ⳱where s is time delay, P is the intensity of the white noise
not depending on frequency.
Orthogonality of terms in the Wiener series with respect
to the Gaussian white noise provides essential advantages indescribing nonlinear systems: (1) effective methods can beused in estimating kernels, based on determination of crosscorrelation functions; and (2) limited numbers of these terms
can represent the best approximations for a real system fromthe viewpoint of the minimal mean square error. By analogy
with an ordinary impulse characteristic h(t), kernel series{hm} can be treated as a generalized, composed impulse
characteristic of a system. The first-order kernel determines
the linear part of the system response, whereas higher-orderkernels describe the interactions between the values of the
input signal in the past with respect to their influence on the
response at present (Marmarelis and Marmarelis, 1978).
In this study, the kernels of the system were found by
the cross-correlation method described in (Marmarelis andMarmarelis, 1978). The zero-order kernel is equal to the average of the response of the system to the Gaussian white
noise:

冮 h (s)x(tⳮs)ds,(5)1y1(t) ⳱ y(t) ⳮ G1[h1, x(t)] ⳮ G0[h0, x(t)],

(11)

0

and the third-order kernel isG2[h2; x(t)] ⳱冮 冮 h (s ,s )x(tⳮs )x(tⳮs )ds ds
2

1

2

1

0

2

1

h3(s1, s2, s3) ⳱(1/6P3)E[y2(t)x(t ⳮ s1)x(t ⳮ s2)x(t ⳮ s3)],

2



ⳮP

冮 h (s ,s )ds ,
2

1

1

(6)

1

(12)

where

0

y2(t) ⳱ y(t) ⳮ G2[h2, x(t)]ⳮ G1[h1, x(t)] ⳮ G0[h1, x(t)]G3[h3; x(t)] ⳱

(13)



冮 冮 冮 h (s ,s ,s )x(tⳮs )x(tⳮs )x(tⳮs )ds ds ds ⳮ
3

1

2

3

1

2

3

1

2

3

0



ⳮ3P

冮 冮 h (s ,s ,s )x(tⳮs )ds ds ,
3

0

1

2

2

1

1

2

(7)

(Marmarelis and Marmarelis, 1978).
The mean square error of approximation is a quantitative measure of the agreement between the model and thereal system response. The constant component of the response corresponding to the zero-order kernel {h0} gives azero-approximation error. The mean square deviation of the385

Nonlinear Seismic Effects in Soils: Numerical Simulation and Study

response from this line is calculated and normalized so that
it is 100. Then the mean square deviation is estimated fromthe response from the model {h0, h1}, which uses the firstorder kernel, and is normalized by the same value as waspreviously normalized for the zero-approximation error. The
deviations then are calculated and normalized by the samevalue of the real system response from responses of the nonlinear models predicted by the kernels {h0, h1, h2} and {h0,h1, h2, h3} (Marmarelis and Marmarelis, 1978). Performing
these calculations, we estimate the degree to which thesecond-order {h2} and the third-order {h3} kernels contribute to the nonlinear component of the system response, andwe draw conclusions regarding the types and quantitative
characteristics of the system nonlinearity. If the secondorder kernel {h2} gives the largest contribution, the system
possesses mostly quadratic nonlinearity; if the third-orderkernel {h3} contributes more, the system is cubic-nonlinear,etc. The system is virtually linear, if nonlinear components
in the response are negligibly small.
In spectral terms, the response of a nonlinear system canbe represented by the sequence of spectra: the power spectrum, bispectrum, trispectrum, and other higher-order spectra(Marmarelis and Marmarelis, 1978; Nikias and Raghuveer,1987). The Nth-order spectrum CN(x1, x2, . . . , xNⳮ1) of
the real discrete stationary process x(k) is defined as theFourier transform of its Nth-order cumulant sequencecN(s1, s2, . . . sNⳮ1), i.e.,c2(s1) ⳱ m2(s1) ⳱ E[x(k)x(kⳭs1)],(18)

c3(s1,s2) ⳱ m3(s1, s2) ⳱ E[x(k)x(kⳭs1)x(kⳭs2)],(19)(where E[ ] means averaging), this is not true for the fourthorder (and higher) statistics. To generate the fourth-order(and higher) cumulant sequence, we need to know the moment and autocorrelation sequences.
Analysis of higher-order spectra of signals passing
through a nonlinear system also allows determination of thetypes and quantitative parameters of the nonlinearity of thesystem (Nikias and Raghuveer, 1987).
In this work, normalized higher-order spectra, that is,higher-order coherences are estimated. One-dimensional diagonal (x1 ⳱ x2 ⳱ . . . ⳱ xNⳮ1 ⳱ x) fourth-, fifth-, andsixth-order coherences are constructed by analogy with thebicoherence r 23(x1, x2) (Haubrich, 1965):
r 32(x1, x2) ⳱|C3(x1, x2)|2.
C2(x1)C2(x2)C2(x1 Ⳮ x2)

兺 . . . s 兺⳱ⳮ cN(s1, s2, . . . , sNⳮ1) .
s ⳱ⳮ
1

(14)

Nⳮ1

• exp{ⳮj(x1 s1 Ⳮ . . . Ⳮ xNⳮ1 sNⳮ1)}Ⳮ兺 c2(s1) exp{ⳮj(x1s1)},s ⳱ⳮ

r 42(x) ⳱r 52(x) ⳱bispectrum (N ⳱ 3):

兺

(16)
Ⳮ

兺 c3(s1, s2) exp{ⳮj(x1s1 Ⳮ x2s2)},m⳱ⳮ n⳱ⳮ

trispectrum (N ⳱ 4):
C4(x1, x2, x3) ⳱







兺 兺 兺 c4(s1, s2, s3) •s ⳱ⳮ s ⳱ⳮ s ⳱ⳮ
1

2

3|C5(x)|2,
C 42(x)C2(4x)(22)

and the sixth-order coherence is calculated by the equationr 62(x) ⳱

Ⳮ

(21)

(15)

1

C3(x1, x2) ⳱|C4(x)|2,3
C2(x)C2(3x)the fifth-order coherence is calculated by the equation(Nikias and Raghuveer, 1987). The power spectrum, bispectrum, and trispectrum are special cases of the Nth-order spectrum, i.e., power spectrum (N ⳱ 2):
C2(x1) ⳱(20)

The fourth-order coherence (tricoherence) is calculated by
the equationCN(x1, x2, . . . , xNⳮ1)
⳱

Here c2(s1), c3(s1, s2), and c4(s1, s2, s3) are sequences of thesecond-order, third-order, and fourth-order cumulants of theprocess x(k), respectively. While the second-order and thirdorder cumulants and the second-order and third-order moments are identical:(17)• exp{ⳮj(x1s1 Ⳮ x2s2 Ⳮ x3s3)}.

|C6(x)|2.
5
C 2(x)C2(5x)(23)

Higher-order coherences take values from 0 to 1 andhave a simple physical meaning, which can be explained by
the following example. Let a signal at the input of a systemwith quadratic (or cubic) nonlinearity contain, among others,
oscillations at frequencies x1 and x2. Then the output signal
will contain oscillations at frequencies x1, x2, as well as at
combination frequencies 2x1, 2x2, x1 Ⳮ x2, and x1 ⳮ x2
(in the case of cubic nonlinearity, 3x1, 3x2, x1 Ⳮ x1 Ⳮ x2,
x1 Ⳮ x2 Ⳮ x2, x1 Ⳮ x1 ⳮ x2, etc.), the phases of thesesignals are connected by certain relations called quadratic(cubic) phase coupling. For two given frequencies, x1 andx2, bicoherence r 23(x1, x2) is equal to the part of the power

386
spectral density at frequency x1 Ⳮ x2, which is coupled withharmonics at frequencies x1 and x2. For three given frequencies x1, x2 and x3, respectively, tricoherencer 24(x1, x2, x3) is the part of the power spectral density at
frequency x1 Ⳮ x2 Ⳮ x3, coupled with harmonics x1,
x2, and x3. If the values of higher-order coherencer 2N(x1, x2, . . . , xNⳮ1) are close to 1, a high phase coherence takes place. Therefore, the combination-frequency harmonic x1 Ⳮ x2 Ⳮ . . . Ⳮ xNⳮ1 is the result of the interactionof the main harmonics x1, x2, . . . , xNⳮ1. On the other side,
if the values are close to 0, the harmonics x1, x2, . . . , xNⳮ1
and x1 Ⳮ x2 Ⳮ . . . Ⳮ xNⳮ1 are independent.
In the considered example, the generation of harmonics
at frequencies 2x1, 2x2, x1 Ⳮ x2, and x1 ⳮ x2 proceeds
simultaneously as a common process, and the knowledge of
one of these processes (for example, generation of harmonicx1 Ⳮ x2) is sufficient for getting a representation of thenonlinearity of the system and of other processes. Respectively, higher-order spectra possess symmetry properties,
and their different symmetric regions correspond to thesedifferent harmonics. Knowledge of a higher-order spectrum
in one of its symmetric regions (for example, knowledge ofthe bispectrum in the triangular region x2 ⱖ 0, x1 ⱖ x2,
x1 Ⳮ x2 ⱕ p) is enough for its complete description (Nikias
and Raghuveer, 1987), as well as for getting a representationof nonlinear interactions in the medium. A special case is
the generation of harmonics at frequencies 2x1 and 2x2,
which correspond to diagonal values of higher-order spectra.
Without essential loss of information, we can use the more
descriptive one-dimensional higher-order spectra instead ofclumsy multidimensional higher-order spectra to get a representation of the nonlinearity of the system.
In this work, conventional direct methods for higherorder spectral estimation were used based on fast Fouriertransform (FFT), which gave nonbiased and consistent estimates. To estimate higher-order spectra for the sequencex(k), the sequence is divided into p equal nonoverlapping
intervals, over which averaging of estimates is performed.
Estimates become smooth for sufficiently large p (Nikias
and Raghuveer, 1987). If x(k) is a stationary Gaussian process, for large p, bicoherence values r 23(x1, x2) are approximately v2 distributed with 2 degrees of freedom and the
expected mean value of 1/p for nondiagonal terms (Haubrich, 1965). It can be easily shown that bicoherence, tricoherence, fifth-order coherence, and sixth-order coherence values are approximately v2 distributed with 2 degrees of
freedom, and the expected mean value of (N ⳮ 1)!/p for
diagonal terms, where N is the order of the coherence, andthe multiplier (N ⳮ 1)! appear because of the symmetry
properties of higher-order spectra. For a Gaussian process,which passed through a nonlinear system, mean values ofhigher-order coherence are larger than (N ⳮ 1)!/p. In thiscase, high values of bicoherence testify to the quadratic nonlinearity of the system and give its quantitative estimates.
High values of tricoherence (or the fifth- and sixth-ordercoherences) testify to the cubic (or fifth- and sixth-order)

O. V. Pavlenko

nonlinearity of the system and give its quantitative estimates.
Thus, we determine the types and quantitative parameters of
the system nonlinearity.
In engineering seismology, the propagation of a vertically incident seismic wave in horizontally layered sediments is also a problem of the transformation of the seismic
wave by a nonlinear system. Here, the input signal is theincident seismic wave, and the output signal is the movementon the surface. In this work, numerical simulation was performed, and the movement on the surface was calculated forthe input signals of three types: the Gaussian white noise,monochromatic signals, and sums of two monochromatic
signals. First, the system was tested by the Gaussian white
noise, and the Wiener kernels were calculated by the methodof cross-correlation functions (Marmarelis and Marmarelis,1978); linear parts of the ground response and nonlinear corrections were determined. Then monochromatic signals at
frequencies of 2, 5, 8, 11, 14, and 17 Hz and pairs of thesesignals were used as input signals, and their nonlinear distortions in the medium were studied. Input signals of various
intensities were used to establish regularities in nonlinear
seismic processes in soils.
To calculate oscillations on the surface, an algorithm
was used elaborated by Joyner and Chen (implicit method,
Joyner and Chen, 1975), in which the rheological nonlinear
elastic-plastic Iwan model was applied, which is similar tothat obtained in laboratory experiments by Hardin andDrnevich (1972). Hysteresis loops shown in Figure 1a correspond to this model. The accuracy of the modeling depends upon the number n of elements used. Calculations
were performed for n ⳱ 50 as in the original article (Joyner
and Chen, 1975). This algorithm implies the stress-strain
relationship to be a transfer function of the medium, which
was very satisfactory for the purposes of this study.
The soil profile used in the numerical simulation was
the same as in the article by Joyner and Chen (1975): ahorizontally layered sedimentary column bounded above by
the free surface and below by a semi-infinite elastic medium
representing the bedrock. The same normalized curves
shown in Figure 1a were used for all of the soil layers; differences in soil behavior from one layer to another resulted
from differences in the values of shear stress in failure smax
and the low-strain shear modulus Gmax assigned to different
layers. The soil profile was a 200-m section of firm alluvium
with a density of 2.05 g/cm3. Shear stress in failure, smax,
and the S-wave velocity progressively increased with depth.
Parameters of the underlying half-space were the same as in
the article by Joyner and Chen: the substratum was assigned
a shear velocity of 2.0 km/sec, and a density of 2.6 g/cm3;
part of the energy of the incident wave was radiated backinto the underlying medium.
The choice of the nonlinear model depended on a few
factors. On the one hand, this model captured common elements of nonlinear soil behavior that made it useful for abaseline evaluation of site response. On the other hand, it
was simple, and therefore it was advantageous for studying

Nonlinear Seismic Effects in Soils: Numerical Simulation and Study

new phenomena. In the study presented here, the complexityof the model could be a nuisance, and it is shown below thatthe obtained conclusions are of a rather common character.
The choice of the nonlinear model influences quantitative
characteristics of nonlinearity, but cannot add anything qualitatively new to the results. The same can be said about thesite geometry. Only quantitative estimates are influenced by
the site geometry. According to the theory of wave propagation in a nonlinear medium, nonlinear distortions are accumulated, therefore, the thicker the sedimentary layer, themore pronounced the manifestations of nonlinearity (Zarembo and Krasil’nikov, 1966).

Results and Discussion
Figure 2a, b, and c represents examples of the responseof the studied system of soil layers to the Gaussian noise ina wide frequency range (0.3–170 Hz), obtained in numericalsimulation and predicted by models. Estimates of the zero-,first-, second-, and third-order kernels of the system werebased on input and output signals of 65,000 points of duration (190 sec).
Zero-order kernels h0 are constant components of theseismic wave field at the surface, namely, the quasi-static
deformations of the surface. They are small, but nonzero;for different time intervals of the response they take positiveor negative values, which are larger for higher-intensity input signals. The limits of their variations are approximately
proportional to the mean square amplitudes V2(t) of oscillations on the surface, which is in agreement with the theoryof nonlinear wave interactions (Zarembo and Krasil’nikov,
1966), and they are presented in Table 1. This means thathorizontal displacements of points on the surface in positive
and negative directions are different, that is, subsurface layers slowly shift to one or the other side in their oscillations.
This effect can be interpreted as a result of accumulation ofresidual shift deformations on the surface. Asymmetry,skewness of oscillations, is concerned with even nonlinearities, which were fairly small in this case, considering thesmall values of h0.
First-order kernels h1(s) were estimated in 300 points,
two-dimensional second-order kernels h2(s1, s2) were estimated in 7225 points, and three-dimensional third-order kernels h3(s1, s2, s3) were estimated in 125,000 points. Basedon these kernels, the response of the linear model (Fig. 2a–c,
trace 3) and nonlinear corrections accounting for quadratic(Fig. 2a–c, trace 5) and cubic (Fig. 2a–c, trace 7) nonlinearities were constructed according to equations (3) to (7).
Trace 4 (Fig. 2a–c) represents the deviation of the real system response from the response of the linear model constructed by the zero- and first-order kernels {h0, h1}. Trace
6 (Fig. 2a–c) is the deviation of the system response fromthe response of the nonlinear model predicted by kernels {h0,h1, h2}. Trace 8 (Fig. 2a–c) is the deviation of the system
response from the response of the nonlinear model predictedby kernels {h0, h1, h2, h3}. For the studied system, consid-

387
ering quadratic and cubic nonlinearities decreases the meansquare error of approximation as is shown in Table 1.
Progressive decrease in the mean square error showedthat the system is substantially nonlinear, and the secondand third-order kernels considerably improved the system’s
description. Accounting for quadratic nonlinearity decreased
the mean square error only by 2%, whereas accounting for
cubic nonlinearity decreased it by 6–10%, that is, the system
possessed a weak quadratic nonlinearity and a substantial
cubic nonlinearity. Note that the intensity of the quadratic
and cubic nonlinear components increases with the increasing intensity of the input signal.
A fairly large residual remained after subtracting thelinear component and the corrections due to the quadratic
and cubic nonlinearities from the calculated system response. The residual can be related to the influence of thehigher-order (fourth-order, fifth-order, etc.) nonlinearities, aswell as to inaccuracies in estimating kernels. These inaccuracies are caused by various reasons, such as the finite duration and the boundedness of the spectral band of the Gaussian noise, truncation of the tails of normal distribution ofnoise amplitudes, and inaccuracies in calculations (the detailed error analysis and correction methods are given in
Marmarelis and Marmarelis, 1978).
Analyzing changes in shapes and spectra of seismic
waves (Fig. 3a–c), three cases are considered: (a) the input
signal is monochromatic, (b) it is a sum of two monochromatic signals, and (c) it is the Gaussian noise in a wide
frequency range (0.3–170 Hz). The first row of seismograms(1) is the input signal, the second row shows the oscillations
on the surface for four different intensities of input signals
corresponding to those listed in Table 1, and the bottom row(8) shows the respective stress-strain curves. Increase inamplitudes of input signals shifts working intervals of thecurves to nonlinear domains, which induces nonlinear distortions: changes in amplification and spectral compositionof signals and acquisition of high-frequency and lowfrequency combination harmonics, namely, nonlinear effects(Zarembo and Krasil’nikov, 1966).
For a nonlinear system, the spectrum of the output signal always differs from the spectrum of the input signal;
usually, it is more complicated. Changes in the spectral composition of signals are shown in row (3) in Figure 3a–c. Input
signal spectra are shown with thin lines, and corresponding
spectra of output signals on the surface are shown with thicklines. Given that the input signal is monochromatic, its nonlinear distortions are due to generation and growth of thethird, fifth, seventh, and so on, harmonics of the main frequency (Fig. 3a: higher harmonics at frequencies of 15, 25,
35, 45, and 55 Hz); therefore, we can conclude that the system possessed odd types of nonlinearity.
The richer the spectrum of the input signal is, the more
complicated the spectrum of the response of a nonlinear system is, which is a result of wave interactions that are morediverse. The superposition principle is not satisfied in a nonlinear system; the system response to a sum of two signals

388

O. V. Pavlenko

Figure 2.

The response of the system of soil layers to theGaussian white noise of various intensities: (a) ⳮ冪V 2(t) ⳱
0.001 m/sec, (b) ⳮ冪V 2(t) ⳱ 0.005 m/sec, (c) ⳮ冪V 2(t)
⳱ 0.02 m/sec. 1 is the input signal; 2 is the response of thesystem; 3 is the response predicted by a linear model; 4 isthe difference between the response of the system and the
response predicted by a linear model; 5 is the nonlinear correction due to quadratic nonlinearity predicted by kernel{h2}; 6 is the difference between the response of the systemand the response predicted by the model {h0, h1, h2}; 7 isthe nonlinear correction due to cubic nonlinearity predictedby kernel {h3}; 8 is the difference between the response of
the system and the response predicted by the model {h0, h1,
h2, h3}.

Nonlinear Seismic Effects in Soils: Numerical Simulation and Study

Table 1
Limiting Values of h0 and Mean Square Errors of Approximation
Intensity 冪V2(t) (m/sec)

Mean Square Errors for Models

Input
SignalOutputSignalLimiting
h0 (m/sec){h0}{h0, h1}{h0, h1, h2}{h0, h1, h2, h3}

0.001
0.005
0.02
0.1

0.003
0.01
0.03
0.1

10ⳮ6
10ⳮ5
10ⳮ4
10ⳮ3

100
100
100
100

32.4
33.8
41.4
56.7

30.4
32.5
38.9
54.0

25.6
26.7
33.2
43.3

differs from a sum of the responses to each of these signals.
If an input signal is a sum of two harmonics, for example,
at frequencies of 5 and 11 Hz (Fig. 3b), its propagation in anonlinear medium is accompanied by the generation ofhigher third, fifth, seventh, and so on, harmonics of each ofthese frequencies: 15, 25, 35, 45, 33, 55 Hz, etc., as well as
by the generation of combination-frequency harmonics of
21 Hz (5 Ⳮ 5 Ⳮ 11 Hz), 27 Hz (5 Ⳮ 11 Ⳮ 11 Hz), 17 Hz(11 Ⳮ 11 ⳮ 5 Hz), 37 Hz (5 Ⳮ 5 Ⳮ 5 Ⳮ 11 Ⳮ 11 Hz),
43 Hz (5 Ⳮ 5 Ⳮ 11 Ⳮ 11 Ⳮ 11 Hz), 39 Hz (11 Ⳮ 11 Ⳮ
11 Ⳮ 11 ⳮ 5 Hz), 49 Hz (11 Ⳮ 11 Ⳮ 11 Ⳮ 11 Ⳮ 5 Hz),and many other combination-frequency harmonics of thefifth, seventh, and so on, orders as a result of higher-orderinteractions. At the same time, combination-frequency harmonics, evidently indicating even types of nonlinearity, areabsent in the spectra of output signals.
Generation of combination-frequency harmonics leads
to redistribution of energy between spectral components. For
rather intense input seismic signals, the power spectral
density of oscillations on the surface takes the form ofE(f )  fⳮk (where f is frequency and k depends on the properties of the medium), regardless of the spectral compositionof an input signal (Fig. 3a–c, the last columns). A similar
dependence was derived theoretically in nonlinear acoustics.
Nonlinear wave effects in the acoustic field cause redistribution of energy over the spectrum. In the case of the interaction of a large number of waves, the power spectral densityof the resulting field is characterized by the dependence ofE(f )  f ⳮ2, provided that the medium is nondispersive anddoes not contain sources and absorbers of energy (Kadomtsev and Karpman, 1971). Dispersion and absorption hinder
the development of nonlinear processes, that is, dispersion
impedes the transition of energy of the main-frequency harmonic to its higher harmonics, whereas absorption leads toa quicker attenuation of higher harmonics. Therefore, coefficient k depends on dispersion, absorption, and nonlinear
properties of the medium.
Therefore, the results of the simulation show that generation of combination-frequency harmonics is observed
even for small amplitudes of input signals, whereas othernonlinear effects, such as decrease in amplification and decrease in the shear modulus, become noticeable only forrather large input signals, when changes in the spectra of
propagating signals achieve their maximum (Fig. 3a–c, rows
2, 3, and 8). Fig. 3a–c, row 8 represents the working areas389
of stress-strain curves for various input signals. These figures
show that decrease in amplification and decrease in the shear
modulus are not noticeable for small amplitudes of input
signals.
Based on the results of the simulation, we concluded
that increase in high-frequency components is expected incases of weak nonlinearity (i.e., weak input signals and/or
thin sedimentary deposit) but is not expected in cases ofstrong nonlinearity (i.e., strong input signals and/or thicklayer of sediments). In the latter case, input signals with any
arbitrary spectra will be transformed into output signals on
the surface with spectra of the type of E(f )  fⳮk, where kdepends on the properties of the medium. A demonstrative
example of such spectral changes are the acceleration records of the 1995 Kobe earthquake at Port Island downhole
group, where high-frequency spectral components are decreasing and low-frequency components are increasing inaccelerogramms from the depth to the surface (Aguirre andIrikura, 1997).
Power spectra and higher-order coherences of input andoutput signals characterize quantitatively nonlinear distortions of signals in the medium. Figure 3a–c represents bicoherence r 23(x) (row 4), tricoherence r 24(x) (row 5), fifthorder coherence r 25(x) (row 6), and sixth-order coherencer 26(x) (row 7). Thin lines correspond to the higher-ordercharacteristics of input signals, and thick lines correspond to
those for surface oscillations.
Figure 3a, b, and c (rows 4, 5, 6, and 7) show that increase in the intensity of input signals leads to increase intricoherence and the sixth-order coherence of signals on thesurface. If an input signal is monochromatic or a sum of
monochromatic signals (Fig. 3a, b), increase in the fourthand sixth-order coherences at the main frequencies of the
input signal (i.e., at 5 Hz and 11 Hz) indicates the generationof their third and fifth higher harmonics. For input signal
amplitudes that are rather high, the fourth- and sixth-ordercoherences are close to 1 at the main frequencies of input
signals (5 Hz—Fig. 3a; 5 Hz and 11 Hz—Fig. 3b). Thismeans that an entire phase coupling occurs between the mainfrequency and its third and fifth harmonics.
The results show that the generation of higher harmonics is more effective for higher main frequencies (i.e., more
effective for frequency 11 Hz than for 5 Hz, Fig. 3b). As
was mentioned previously, propagation of monochromatic
signals at frequencies of 2, 5, 8, 11, 14, and 17 Hz wasstudied, and although working intervals of the stress-strain
diagrams were similar, the third and fifth harmonics weregenerated most effectively for the input signal at frequency
17 Hz. This is in agreement with the theory of nonlinear
wave interactions: amplitudes of combination-frequency
harmonics are proportional to their wave numbers, that is,higher-frequency combination harmonics are generated
more effectively than lower-frequency ones (Zarembo andKrasil’nikov, 1966). The conclusion of a more effective generation of higher harmonics for higher frequencies of input
signals, as well as the alternative sign of quasi-static shift

390

O. V. Pavlenko

Figure 3.

Caption on page 393.

Nonlinear Seismic Effects in Soils: Numerical Simulation and Study

Figure 3.

(Continued).

391

392

O. V. Pavlenko

Figure 3.

(Continued).

393

Nonlinear Seismic Effects in Soils: Numerical Simulation and Study

deformations on the surface, agree with the results describedin Zvolinskii (1982), in which the propagation of a sinusoidal wave in an elastic or plastic layer was calculated. Positive
and negative quasi-static deformations of the surface were
measured in field experiments with a vibrator seismic source
at distances of hundreds of meters from the vibrator (Nikolaev et al., 1995).
If amplitudes of an input sinusoidal signal are rather
high, its propagation in a nonlinear medium is accompaniedby large distortions; the number of combination-frequency
harmonics is high, and the signal becomes noiselike (Fig.
3a,b, last columns). The power spectral density takes theform of E(f )  fⳮk, and the mean values of higher-ordercoherences approach the value of (N ⳮ 1)!/p, which is characteristic for noiselike signals, and exceed it. If the intensityof the input signal is high, higher-order coherences increasein their high-frequency domains (Fig. 3b,c, last columns)
and become greater than (N ⳮ 1)!/p, indicating an essential
content of combination-frequency harmonics in the response. Table 2 represents the mean values of higher-ordercoherences (averaged over the spectral range) of oscillations
on the surface for the input signal–Gaussian noise of various
intensities. The lower row shows corresponding values(N ⳮ 1)!/p (p ⳱ 55 in the calculations).
We can see from these data and from Figure 3c that if
the intensity of the input Gaussian noise is high, all higherorder coherences exceed values of (N ⳮ 1)!/p. However,calculated real and imaginary components of the higherorder coherences (Fig. 4) show that increase in the third- andfifth-order coherences is caused by the increase in scatteringof their real components about a zero mean, whereas the realcomponents of the fourth- and sixth-order coherences actually increase in the high-frequency range. Corresponding
imaginary components vary about zero with zero means forany amplitudes of input and output signals. Therefore, increase in the fourth-order and sixth-order coherences indicates the cubic and the fifth-order nonlinearities of thesystem.
Thus, the numerical simulation shows that subsurface
soils with a typical stress-strain relationship (Fig. 1a) arecharacterized mostly by odd types of nonlinearity, whereas
even nonlinearities are weak. This is quite natural because
the stress-strain curves describing the conversion of input
signals into the response are described by odd functions withweak even components; in this case, even-order spectra of

Table 2
Mean Values of Higher-Order Coherences of Oscillations
on the Surface
Intensity of the
Input Signal 冪V2(t) (m/sec)

r 23(x)

r 24(x)

r 25(x)

r 26(x)

0.001
0.005
0.02
0.1
(N ⳮ 1)!/p

0.037
0.04
0.041
0.066
0.036

0.105
0.107
0.119
0.35
0.11

0.38
0.43
0.48
0.99
0.43

1.8
1.72
2.12
4.49
2.1

output signals were close to zero (Nikias and Raghuveer,1987). For higher amplitudes of input signals, working intervals of the stress-strain curves shift to nonlinear areas,
where even components gain in magnitudes. In these cases,
first of all, quasi-static deformations of the surface appear asa summed effect of all even-order nonlinearities. Even-ordernonlinearities are characteristic for media possessing stressstrain dependencies, which are described by functions withsignificant even components. For example, such behavior isobserved in liquefied soils and in cases in which the stressstrain curves contain parts with large slopes (close to horizontal). In such media, additionally to quasi-static deformations, generation of even-order higher harmonics can be
expected.
Quantitative manifestations of nonlinearity, that is, amplitudes of combination-frequency harmonics, are definedby the stress-strain relationships. As is known, real stressstrain dependencies of soils obtained in experiments are very
diverse and dependent on the granulometric composition ofa soil, its humidity, and so on. Stress-strain curves for watersaturated soils substantially differ from the stress-strain relationship shown in Figure 1a: the initial convex-up part isfollowed by a deviation to the stress axis (concave-up part)(Vasil’ev et al., 1977). The character of the seismic wave
field and the effects of nonlinearity are different in such
soils, and we can expect shock waves beyond the limit ofelasticity (Zvolinskii, 1982). Stress-strain dependencies represented by the only convex-up parts are sometimes calledsoft diagrams (this type is characteristic for soft soils likeloess loams and sands); whereas so-called “hard” dependencies composed of concave-up parts are characteristic forwater-saturated sands and clays. Respectively, different indications of nonlinear behavior are known in sands and clays:Figure 3.

Changes in shapes, spectra, and higher-order coherences of seismic signals propagating in a nonlinear medium (three cases are considered, with different
intensities of input signals): (a) a monochromatic signal, (b) a sum of two monochromatic signals, (c) the Gaussian noise in a wide frequency range. Row 1, input signals;
row 2, the response (oscillations on the surface); row 3, power spectra of the input
signals (thin lines) and spectra of the responses (thick lines); rows 4–7, one-dimensional
bicoherence r 23(x), tricoherence r 24(x) , fifth-order r 25(x) , and sixth-order r 26(x) coherences of the input signals (thin lines) and the responses (thick lines); row 8, corresponding stress-strain diagrams (they depend on depth; maximum values are shown).

394

O. V. Pavlenko

Figure 4.

Real (rows 1–4) and imaginary (rows 5–8) components of higher-ordercoherences of the Gaussian noise in a wide frequency range: rows 1, 5, real and imaginary components of bicoherence; rows 2, 6, real and imaginary components of tricoherence; rows 3, 7, real and imaginary components of the fifth-order coherence; rows 4,
8, real and imaginary components of the sixth-order coherence. Thin lines show characteristics of input signals. Thick lines show that of output signals.

395

Nonlinear Seismic Effects in Soils: Numerical Simulation and Study

nonlinear effects are more often observed in sandy sites,
whereas nonlinear behavior of clay deposits with high plasticity index appears at higher deformations (Bard, 1995).
To estimate variations in quantitative characteristics ofnonlinearity for different types of soils, the calculations wererepeated for cases of a more soft stress-strain relationshipthan that obtained by Hardin and Drnevich (1972), and fora more hard type of the relationship as well. The results werequalitatively the same: the higher harmonics of odd types
were generated. For weak input signals, variations in theamplitudes of higher harmonics can achieve 200–300% depending on the model. For high amplitudes of input signals,
the exponent k differs for different models.

Conclusions
Numerical simulation was performed of the propagationof vertically incident seismic waves in a system of horizontal
soil layers. The method of nonlinear system identification
was applied to analyze nonlinear properties of soils. Soils
possessing a typical stress-strain relationship obtained in laboratory experiments by Hardin and Drnevich (1972) werestudied.
The results of the numerical simulation show that nonlinear wave processes develop in soils in the following order:
the third, fifth, seventh, and so on, higher harmonics of themain frequencies of the input signal are generated, and theinteraction between them leads to generation of new combination-frequency harmonics. As a result, the number ofspectral components increases, and the spectrum takes theform of E(f )  fⳮk, in which k depends on dispersion, absorption, and nonlinear properties of the medium. Generation of higher harmonics was found to be more effective forhigher-frequency input signals.
Therefore, in cases of weak nonlinearity (i.e., weak input signals and/or thin sedimentary layer), an increase inhigh-frequency components is expected, whereas in cases ofstrong nonlinearity (i.e., strong input signals and/or thicksedimentary layer), input signals with arbitrary spectra aretransformed into output signals on the surface with spectraof the type of E(f )  fⳮk, where k depends on the propertiesof the medium. A weak nonlinearity corresponds to the casein which manifestations of nonlinearity such as degradationof shear moduli and decrease in signal amplification are notnoticeable, whereas, strong nonlinearity corresponds to thecase in which decrease in amplification and in shear moduli
can be detected.
In real media, dispersion, absorption, and seismic noisehinder the development of nonlinear processes. As a result,
the third- and fifth-order nonlinearities should be of mostimportance in real media, whereas more higher-order nonlinearities, such as the seventh-, ninth-order are apparently
negligible. The influence of even-order nonlinearities andthe possibility of static deformations seem to be small inusual cases, however, they become significant in cases, whena stress-strain relationship of a soil is described by functions

with noticeable even components (for example, a stressstrain dependence containing parts with slopes close to horizontal).
As a whole, nonlinear seismic processes in soils are determined by their stress-strain dependencies, by the intensity
and spectral composition of the incident seismic wave, andby dispersion and absorption properties of the medium.

Acknowledgments
I wish to thank my colleagues from Radiophysical Research Institutein Nizhnii Novgorod, Vladimir V. Gushchin for his valuable comments anddiscussions, and Vladimir A. Shemagin for helpful recommendations concerning higher-order spectral estimation. I also thank Aleksei V. Nikolaev
for his careful supervision to this research. Significant improvement to thearticle resulted from the comments of Arthur F. McGarr, Jonathan P. Stewart, and an anonymous reviewer. I am particularly grateful to Professor
Kojiro Irikura for his encouragement and support of this work. This study
was partially supported by Grants-in-Aid for Scientific Research, Number
11792026 and Number 00099 (JSPS), from the Ministry of Education, Science, Sports, and Culture of Japan.

References
Aguirre, J., and K. Irikura (1997). Nonlinearity, liquefaction, and velocity
variation of soft soil layers in Port Island, Kobe, during the 1995
Hyogo-ken Nanbu earthquake, Bull. Seism. Soc. Am. 87, no. 5, 1244–
1258.
Aki, K., and K. Irikura (1991). Characterization and mapping of earthquake
shaking for seismic zonation, in Proc 4th Int. Conf. on Seismic Zonation, Stanford, California, Vol. 1, 25–29, August, 61–110, EERI,Oakland, California.
Aleshin, A. S., V. V. Gushchin, M. M. Krekov, A. V. Nikolaev, and G. M.
Shalashov (1981). Experimental studies of nonlinear interaction ofseismic surface waves, Dokl. Akad. Nauk SSSR, 260, no. 3, 574–575.
Bard, P.-Y. (1995). Effects of surface geology on ground motion: recent
results and remaining issues, in Proc. of the 10th European Conf. on
Earthquake Engineering, G. Duma (Editor), Vol. 1, Vienna, August–
September 1994.
Bard, P.-Y., and K. Pitilakis (1995). Seismic zonation and ground motion
interface, in Proc. of the 5th International Conference on Seismic
Zonation, Vol. III, Nice, France, 17–19 October, Ouest Editions,
Presses Academiques, Nantes, France, 2127–2152.
Beresnev, I. A., Wen K.-L., and Y. T. Yeh (1995). Nonlinear soil amplification: its corroboration in Taiwan, Bull. Seism. Soc. Am. 85, 496–
515.
Caillot, V., and P.-Y. Bard (1990). Characterizing site effects for earthquake regulations in the French seismicity context: a statisticalanalysis, in Earthquake Soil Movement; Bases and Foundations;
Structure-Soil Interaction; Models of Seismic Forces Influence, Proc.
of the 9th European Conf. on Earthquake Engineering, Vol. 4B, Moscow, 27–36.
Chang, C.-Y., C. M. Mok, and H. T. Tang (1996). Inference of dynamic
shear modulus from Lotung downhole data, J. Geotech. Eng. 122,
no. 8, 657–665.
Darragh, R. B., and A. F. Shakal (1991). The site response of two rock andsoil station pairs to strong and weak ground motion, Bull. Seism. Soc.
Am. 81, 1885–1899.
Dimitriu, P. P. (1990). Preliminary results of vibrator-aided experiments innonlinear seismology conducted at Uetze, F.R.G., Phys. Earth Planet.
Interiors 63, 172–180.
Elgamal, A.-W., M. Zeghal, H. T. Tang, and J. C. Stepp (1995). Lotung
downhole array. I: Evaluation of site dynamic properties, J. Geotech.
Eng. 121, no. 4, 350–362.

396
Groshkov, A. L., R. R. Kalimulin, G. M. Shalashov, and V. A. Shemagin(1990). Nonlinear borehole logging with the method of acoustic wave
modulation by seismic fields, Dokl. Akad. Nauk SSSR 313, no. 1,
63–65.
Gushchin, V. V., and G. M. Shalashov (1981). On the possibility of application of nonlinear seismic effects to problems of vibro-radiation of
the earth, in Study of the Earth by Nonexplosion Seismic Sources,
Nauka, Moscow, 144–155.
Hardin, B. O., and V. P. Drnevich (1972). Shear modulus and damping insoils: design equations and curves, Proc. Am. Soc. Civil Eng. J. Soil
Mech. Found. Div. 98, 667–692.
Haubrich, R. A. (1965). Earth noise, 5 to 500 millicycles per second, J.
Geophys. Res. 70, 1415–1427.
Joyner, W. B., and T. F. Chen (1975). Calculation of nonlinear groundresponse in earthquakes, Bull. Seism. Soc. Am. 65, 1315–1336.
Kadomtsev, B. B., and V. I. Karpman (1971). Nonlinear waves, Uspekhi
Fiz. Nauk 103, no. 2, 27–48.
Kamiyama, M. (1992). Non-linear soil amplification identified empiricallyfrom strong earthquake ground motions, J. Phys. Earth 40, 151–174.
Lund, F. (1983). Interpretation of the precursor to the 1960 Great Chilean
Earthquake as a seismic solitary wave, Pageoph 121, no. 1, 17–26.
Marmarelis, P. Z., and V. Z. Marmarelis (1978). Analysis of Physiological
Systems, The White-Noise Approach, Plenum Press, New York.
Nikias, C. L., and M. R. Raghuveer (1987). Bispectral estimation: a digital
signal processing framework, Proc. IEEE 75, 869–891.
Nikolaev, A. V. (1967). Seismic properties of a loose medium, Izv. Akad.
Nauk SSSR Fiz. Zemli 2, 23–31.
Nikolaev, A. V. (1988). Problems of nonlinear seismology, Phys. Earth
Planet. Interiors 50, 1–7.

O. V. Pavlenko
Nikolaev, A. V., O. V. Pavlenko, and A. P. Yakovlev (1995). Quasi-static
vibration-induced deformations of the Earth’s surface and nonlinear
properties of rocks, Physics of the Solid Earth 30, 1023–1030 (English
translation).
Parker, D. F. (1988). Stratification effects on nonlinear elastic surface
waves, Phys. Earth Planet. Interiors 50, 16–25.
Vasil’ev Yu. I., L. A. Ivanova, and M. N. Shcherbo (1969). Measurementof stress and strain in a soil during the propagation of waves from
explosions, Izv. Akad. Nauk SSSR Fiz. Zemli 1, 21–37.
Vasil’ev Yu. I., A. A. Gvozdev, L. A. Ivanova, L. V. Molotova, V. B.
Fomichev, and, M. N. Shcherbo (1977). Mechanical properties of asoft soil under stresses up to (5–10)•105 Pa and selecting the modelof the soil behavior in strong earthquakes, in Seismic Microzonation,
S. V. Medvedev (Editor), Nauka, Moscow, 121–129.
Zarembo, L. K., and V. A. Krasil’nikov (1966). Introduction to Nonlinear
Acoustics, Nauka, Moscow.
Zeghal, M., A.-W. Elgamal, H. T. Tang, and J. C. Stepp (1995). Lotung
downhole array. II. Evaluation of soil nonlinear properties, J. Geotech. Eng. 121, no. 4, 363–378.
Zvolinskii, N. V. (1982). Wave processes in non-elastic media, Problemsof Engineering Seismology 23, 4–19.
Institute of Physics of the Earth, Russian Academy of Sciences
B. Gruzinskaya 10
Moscow 123995, Russia
olga@synapse.ru
Manuscript received 31 March 2000.