Geophys. J. Int. (2005) 160, 539–553

doi: 10.1111/j.1365-246X.2005.02230.x

Identification of the non-linear behaviour of liquefied
and non-liquefied soils during the 1995 Kobe earthquake
Olga V. Pavlenko1 and Kojiro Irikura2
1 Institute of Physics of the Earth, Russian Academy of Sciences, B. Gruzinskaya 10, Moscow 123995, Russia. E-mail: olga@synapse.ru
2 Disaster Prevention Research Institute, Kyoto University, Gokasho, Uji, Kyoto 611-0011, Japan. E-mail: irikura@egmdpri01.dpri.kyoto-u.ac.jp

SUMMARY
A method based on the non-linear system identification technique is suggested for estimating
the contents of the non-linear components, which are the result of quadratic, cubic and higherorder non-linearities, in the ground response during the strong motion. The method is applied
to data from the near-fault zone of the 1995 Kobe earthquake, and the contents of linear
and non-linear components in the ground response, changing with time during the strong
motion, are estimated for Port Island (PI), SGK and TKS sites, located at 2, 6 and 15 km
from the fault plane, respectively. At PI, the non-linear part of the response increased with
developing liquefaction, it was as high as ∼ 40–60 per cent of the intensity of the response.
At SGK and TKS sites, the non-linear components of the response did not exceed ∼40 and∼13 per cent of the intensity of the response, respectively. Odd-order non-linear components
predominated in the soil response, whereas even-order non-linear components increased andbecame comparable with odd-order ones in liquefied soils and in cases of high intensity ofthe strong motion, when the loading parts of the stress–strain relations of the upper layers
gained noticeable even components. As a whole, the contents of odd-order and even-ordercomponents in the soil response are determined by the shapes of the stress–strain relationsin the upper most non-linear soil layers. At the three sites, changes in spectra of earthquake
signals in subsurface soil layers were substantial as a result of the high non-linearity of thesoil behaviour and spectra of signals on the surface tend to take the form of E( f ) ∼ f −k .
The limiting spectral shape was achieved during the Kobe earthquake at PI and SGK sites.
The proposed methods for processing vertical array records allow understanding of seismic
wave transformations in subsurface soils and are useful for predicting the soil behaviour during
future earthquakes.
Key words: non-linear behaviour of soils, non-linear system identification, strong groundmotion, types of soil non-linearity.

1 I N T RO D U C T I O N
The non-linear behaviour of soils in strong ground motion explains
the well-known fact that on soil sites, records of weak earthquakes
qualitatively differ from records of strong ones. Quantitative characteristics of the soil non-linearity and estimates of the ground response at a given site during the strong motion of any arbitrary
intensity are necessary for practical applications in seismology.
However, because the soil response depends on many factors, suchas the composition and thickness of the soft deposits, their saturation with water, the level of the underground water, and the magnitude and frequency content of an earthquake, the problem is rather
complicated, and for a long time it remained as one of the urgent
problems in seismology, inducing controversies and discussions.
During the 1960s–1970s, after the catastrophic earthquakes in
Niigata and Anchorage, laboratory (Hardin & Drnevich 1972) and
C 2005 RAS

field (Vasil’ev et al. 1969) experiments aimed at studying non-linear
elastic properties of soils began. In these experiments, non-linear
hysteretic stress–strain relations were obtained for various types ofsoils. The questions arose, what corrections should be introduced
into state equations of soils and fractured weathered rock to account for their elastic non-linearity? How large are these corrections for various types of soils and what is the threshold of strains,when they become significant? What is the mutual influence of the
non-linearity of the medium, dispersion and absorption on seismic
wavefields?
To estimate characteristics of elastic non-linearity of geological media, field experiments were performed and non-linear wave
effects in seismic fields were studied, such as, interaction of seismic waves, generation of combination-frequency harmonics and theconstant components of the wavefields, seismic solitary waves, inversion of the seismic wave front, etc. (Aleshin et al. 1981; Gushchin

539

GJI Seismology

Accepted 2004 January 5. Received 2003 September 18; in original form 2002 August 1540

O. V. Pavlenko and K. Irikura& Shalashov 1981; Lund 1983; Shalashov 1984; Groshkov &Shalashov 1986; Dimitriu 1988, 1990; Groshkov et al. 1991, etc.).
In these experiments, effective parameters of non-linearity, definedas the ratios of elastic parameters of the third and second, or the
fourth and second orders (i.e. characteristics of the quadratic andcubic non-linearity), were estimated for subsurface soils. The obtained estimates, N = 102 ÷ 104 (Gushchin & Shalashov 1981;
Shalashov 1984; Groshkov et al. 1991) were substantially higher
than those for crystals and polycrystals: N = 1 ÷ 10 (Zarembo &Krasil’nikov 1966). This fact was explained by the model of the nonlinear deformation of an inhomogeneous elastic–viscous medium.
The model phenomenologically accounts for a large number of nonelastic viscous bonds with various relaxation times, which connect
microblocks and macroblocks (Groshkov & Shalashov 1986). It wastheoretically shown that, in conditions of loading or unloading ofsuch a medium, the most soft bonds experience the largest deformations; as a consequence, effective values of parameters of nonlinearity substantially increase (Groshkov & Shalashov 1986). As aresult, non-linear corrections in state equations of geological materials become comparable with linear terms even at strains as small as∼10−7 –10−4 , according to Gushchin & Shalashov (1981), ∼10−8 ,according to Zimenkov & Nazarov (1993), and ∼10−7 , according
to Van Den Abeele & Johnson (1996).
Observational validations of soil non-linearity were difficult tofind in strong motion records, but the strong motion database grewquickly in the last decades and now seismologists can directly observe the presence of non-linear components in accelerograms as
characteristic waveforms (Archuleta 1998). However, it is not often
clear how much elastic non-linearity of soils influences seismic oscillations on the surface. If only records at the surface are available,
it is difficult to distinguish the influence of topography, inhomogeneities in the crust, source directivity, etc. from the effects of soil
non-linearity. The 1994 Northridge earthquake provided an illustrative example of the fact. Combined results from numerous alluvial
recording sites testify that the decrease in the amplification of themain shock on sediments, compared with the amplification of aftershocks, was caused by a significant non-linearity of sediments (Field
et al. 1997). However, another explanation, though rather difficult,with linear wave propagation models is also possible (O’Connell
1999).
Thus, observations of the strong motion in situ and seismic vertical arrays are necessary for understanding non-linear phenomena insubsurface soils (Kudo 1995; Bard & Pitilakis 1995). The behaviourof soils in situ in strong ground motion is not explored enough, andoften it can not be satisfactorily described by the existing models
and properly understood. This is partially a result of the fact that, instrong ground motion, soils do not only change the parameters of
the propagating seismic waves, but they also change their properties. Phase transitions can occur, when soils are liquefied or saturated
with gas bubbles. As known from physics, such multiphase media
possess strong elastic non-linearity. These processes, occurring insoils in strong ground motion, are important for understanding thesoil behaviour and the necessary experimental data can be obtained
from observations in situ.
Networks of borehole instruments (Kik-Net, Japan; DART, Taipei
basin; Southern California Earthquake Center, Los Angeles basin;
Bay Bridges, San Francisco) are developing worldwide. Accumulation of experimental data on strong ground motion in various soil
conditions could provide the answers to virtually all the questions
that are interesting for seismologists. Where in the soil column does
the response become non-linear? At what level of strain (stress) does
the non-linear response occur? How is the non-linear response man-ifested in the observed ground motion? What properties of the soilthat can be measured in situ are most important in affecting nonlinearity? If a soil column experiences the non-linear response, does
it return to a state with the same linear (or nearly linear) response?
At what time period? (Archuleta & Steidl 2001).
Pavlenko & Irikura (2002b, 2003) have shown that strong motion
records, provided by seismic vertical arrays, allow estimation of
the non-linear hysteretic stress–strain relations in the soil layers
at different depths from the surface down to the location of the
deepest device. Records of the 1995 Kobe earthquake were analysedand vertical distributions of stress–strain relations in soil layers,changing with time during the strong motion, were obtained for Port
Island (PI), SGK and TKS vertical array sites, located at distancesof 2, 6 and 15 km from the fault plane, respectively.
The obtained numerical models of the soil behaviour at the threesites were tested by the Gaussian white noise, and the contents of linear and non-linear components in the soil response were estimatedat the beginning, middle and final periods of the strong motion(Pavlenko & Irikura 2002b). Changes in spectra of the propagating signals were studied and it was found that spectra of signals
on the surface tend to take the form of E ∼ f −n , as a result of
mutual interactions of separate spectral components of propagating
seismic signals. This conclusion is evidently important, because itmeans that in cases of strong non-linearity input signals of any arbitrary spectra are transformed in soil layers into signals of the same
spectrum of the type E ∼ f −n , i.e. we loose information about
the spectral compositions of input seismic signals. Deconvolution
and other procedures aimed at reconstructing parameters of input
signals based on surface records are impossible in cases of strong
non-linearity. Therefore, it would be useful to have a representation,
what strong non-linearity means and how high is the influence ofthe soil non-linearity on the ground response at various distancesfrom the fault plane in various types of soils.
Vertical array records of the 1995 Kobe earthquake at PI, SGK and
TKS sites provide all the necessary data for such a kind of analysis.
This paper describes the detailed study of the soil behaviour at thethree sites during the main shock of the Kobe earthquake and thedetailed analysis of transformations of spectra of seismic signals
in the near-fault areas. Estimates obtained in (Pavlenko & Irikura
2002b) are defined more exactly and completely.
2 M E T H O D A N D D AT A
In system analysis, the non-linear identification of a system implies
determination of linear and non-linear domains of the system response and construction of such a mathematical model of a systemthat its seismic response coincides with the response of the real physical system (Marmarelis & Marmarelis 1978). Soil profiles can be
represented as non-linear systems (because of the non-linear stress–
strain relations), transforming input seismic signals into the groundresponse, i.e. movement on the surface. The most quick and effective
method for the non-linear system identification is testing the studied
system with the Gaussian white noise and calculating the Wiener
kernels (Marmarelis & Marmarelis 1978). If an input is the Gaussian white noise, an output can be represented as the Wiener series(Marmarelis & Marmarelis 1978):y(t) =∞
G m [h m (τl , · · · , τm ); x(t  ), t  ≤ t],

(1)

m=0

where Gm are orthogonal functionals, if x(t) is the Gaussian white
noise with a zero mean, {h m (τ 1 , . . . , τ m )} is a sequence of the
Wiener kernels and τ 1 , . . . , τ m are time delays. The first four WienerC 2005 RAS, GJI, 160, 539–553

Non-linear behaviour during the 1995 Kobe earthquake
functionals are:G 0 [h 0 ; x(t)] = h 0 , ∞G 1 [h 1 ; x(t)] =h 1 (τ )x(t − τ )dτ,

(2)
(3)

0

G 2 [h 2 ; x(t)] =  ∞h 2 (τ1 , τ2 )x(t − τ1 )x(t − τ2 ) dτ1 dτ2 ∞−P
h 2 (τ1 , τ1 ) dτ1 ,
0

(4)

0

G 3 [h 3 ; x(t)] =  ∞

h 3 (τ1 , τ2 , τ3 )x(t − τ1 )0

× x(t − τ2 )x(t − τ3 ) dτ1 dτ2 dτ3
 ∞
− 3P
h 3 (τ1 , τ2 , τ2 )x(t − τ1 ) dτ1 dτ2 ,(5)0where τ is time delay and P is the intensity of the Gaussian white
noise not depending on frequency.
Orthogonality of terms in the Wiener series provides substantial advantages in the non-linear system analysis, because the
Wiener kernels can be easily calculated by the method of the crosscorrelation functions and a limited number of terms in the Wiener
series represent the best approximation for a real system fromthe viewpoint of the minimal mean square error (Marmarelis &Marmarelis 1978). By analogy with an ordinary impulse characteristic h(t), kernel series {h m } can be treated as a generalized, composed impulse characteristic of a non-linear system. The first-orderkernel determines the linear part of the system response, whereas
higher-order kernels describe the quadratic, cubic and higher-ordernon-linear corrections. The non-linear corrections express interactions between the values of the input signal in the past with respect
to their influence on the response at present. Analysing the nonlinear components in the output, we can make a judgement about
the types and quantitative characteristics of the system non-linearity:
if the second Wiener functional G 2 provides the largest contribution
to the system response, the system possesses mostly quadratic nonlinearity; if the third functional G 3 contributes more, the system is
cubic non-linear, etc. (Marmarelis & Marmarelis 1978). Knowledgeof stress–strain relations in soil layers in strong ground motion allows calculation of the propagation of testing signals in the studied
soil profiles and the non-linear identification of the soil behaviour(Pavlenko 2001).
We applied these methods to perform the non-linear identification of the behaviour of soils at PI, SGK and TKS vertical arraysites during the 1995 Kobe earthquake. Fig. 1 shows the locationsof these sites, the major principal axes, epicentres of the main shock
and aftershocks and the soil profiles. The procedure of estimation ofstress–strain relations in soil layers is described in detail in (Pavlenko& Irikura 2003). The soil profiles were divided into groups of layers, according to their composition and saturation with water. For
these groups, physically justified types of stress–strain relations
were assumed. Sets of parametric stress–strain curves were generated and item-by-item examination was applied to find groups of
curves showing the best-fitting approximation to the observed data.
To account for temporal changes in the soil behaviour, records weredivided into intervals of 1.5-s duration. Within each interval, thestress–strain relations were assumed to be stationary and to vary
for different intervals. Calculations were performed successively,interval by interval.
The observed and simulated acceleration time histories of themain shock (horizontal components) and 10 successive time inter
C 2005 RAS, GJI, 160, 539–553

541

vals used in our calculations are shown in Fig. 2. Fig. 3 illustrates estimated hysteretic non-linear stress–strain relations, averaged within
the groups of the upper layers at PI, SGK and TKS sites. The whole
set of stress–strain relations, obtained for all the depths, from the surface down to the location of the deepest device, is given in (Pavlenko& Irikura 2003). Because in our calculations stress–strain relations
are represented in their normalized form as suggested by Joyner& Chen (1975), so that the same relations describe the behaviourof different layers within a group, averaging of stress–strain relations within a group means finding average working intervals of thenormalized stress–strain relations within the group.
Thus, numerical models of the soil behaviour during the 1995
Kobe earthquake were constructed for PI, SGK and TKS sites. Reliability of the obtained representations at the three sites is verified
by the following.
(i) A good agreement between the observed and simulated
records and similarity of the stress–strain relations obtained for two
horizontal components (Figs 2 and 3).
(ii) Physical correctness of the description of the process: progressive liquefaction in the upper layers at PI, reduction and recoveryof the shear moduli in the upper layers at SGK and TKS, and stablebehaviour in deeper parts at the three sites.
(iii) An agreement between our estimates and conclusions ofother authors (Kawase et al. 1995; Sato et al. 1996; Kazama et al.
1998).
The obtained numerical models were used for the non-linear identification of the soil behaviour during the 1995 Kobe earthquake at
the three sites, i.e. for estimation of the contents of non-linear components in the soil response at these sites. In the non-linear identification, before calculating Wiener kernels and functionals, the studied non-linear system is usually tested by monochromatic signals,
in order to reveal the types and roughly estimate some quantitative
characteristics of the system non-linearity.
Testing signals used in our calculations, such as monochromatic
signals and the Gaussian white noise, were approximations to the
input seismic signals during the Kobe earthquake. As seen fromformulae (2)–(5), estimates of the Wiener kernels depend on theintensity of the testing Gaussian white noise. To obtain reliable
estimates of the Wiener kernels, we tried to attain conditions ofpropagation of testing signals close to real seismic waves during
the earthquake. We chose the intensity of testing signals so that theworking intervals of stress–strain relations during the propagationof these signals in the soil layers were close to those ones during the
propagation of signals from the earthquake. However, in a layered
medium, propagation features of seismic signals with different spectral compositions are different and, in general, it is not possible tochoose such intensity of testing signals that the working intervals ofstress–strain relations during the propagation of testing signals coincide with those ones during the earthquake. Therefore, we chose thecontrolling parameters, such as mean-square stresses, mean-squarestrains and mean-square slopes of the stress–strain curves, within
the groups of layers, and in choosing the intensities of testing signals, we tried to adjust these parameters for the groups of layers.
Two groups of layers were considered for SGK and TKS sites andfive groups of layers for PI, as shown in Fig. 1.
In our previous study (Pavlenko & Irikura 2002b), only one parameter was used to adjust the intensity of testing Gaussian white
noise signals, such as an average slope of the stress–strain curves in
the groups of the upper layers. However, numerical examples showthat more exact estimates of Wiener kernels are obtained if we use
more parameters. In addition, we substantially increased the duration

542

O. V. Pavlenko and K. Irikura

Figure 1. Locations of vertical array sites Port Island (PI), SGK and TKS around the Osaka bay, the major principal axes, the epicentres of the main shock
and aftershocks (derived from Sato et al. 1996) and the ground profiles and characteristic types of stress–strain relations.

of testing Gaussian white noise signals compared with our previous
study (Pavlenko & Irikura 2002b), because this also leads to moreaccurate and reliable estimates of Wiener kernels and functionals(Marmarelis & Marmarelis 1978).
3 R E S U LT S A N D D I S C U S S I O N
Fig. 4 presents the results of testing the obtained models of the soil
behaviour (East-West components) during the 1995 Kobe earthquake at the three sites with monochromatic signals. Spectra of
input monochromatic signals are shown in lower rows at each site

and above them are spectra of simulated signals at depths of locations of the recording devices for 10 successive time intervals. As
seen from the figure, signals, propagating up to the surface, acquire
higher harmonics of the main frequencies. The higher harmonics ofthe main frequencies, which are generated during the propagationof the monochromatic signals in the soil layers, indicate the types ofelastic non-linearity of the soil profiles. A large number of the higher
harmonics at PI and SGK sites testifies to a strong non-linearity ofthe soil behaviour at these sites. Whereas at TKS site, few higher harmonics are generated indicating a weaker non-linearity of the soil behaviour at TKS site. Generation of the third, fifth, seventh and otherC 2005 RAS, GJI, 160, 539–553

Non-linear behaviour during the 1995 Kobe earthquake

543

Figure 2. Acceleration time histories of the main shock of the 1995 Kobe earthquake, observed and simulated, at PI (a), SGK (b) and TKS (c) sites. Dashlines indicated selected intervals, for which the results of the non-linear identification are shown in Figs 5(a)–(c).
C 2005 RAS, GJI, 160, 539–553

544

O. V. Pavlenko and K. Irikura

Figure 3. Estimated stress–strain relations in the soil layers during the 1995 Kobe earthquake, averaged within the groups of layers at (a) PI (1: 0–13 m and2: 13–17.5 m), (b) SGK (0–11 m) and (c) TKS (0–14 m), for 10 successive time intervals. Scales in relative units are the same for all time intervals. Stress andstrain are scaled in the manner described by Joyner & Chen (1975): stress is normalized by multiplying by 1/τ max and strain is normalized by multiplying by
G max /τ max .

odd-order harmonics testifies to odd-order non-linearities of the soil
behaviour. Even-order higher harmonics, such as the second, fourth,etc., are generated only in liquefied soils at PI and at SGK site during
the second and third time intervals of the strong motion.
Because non-linear hysteretic stress–strain relations determine
transformations of input seismic signals into the ground response,
generation of higher harmonics obviously depends on the shapes ofthe stress–strain curves. Loading and unloading parts of the hysteretic stress–strain curves can be represented as sums of even andodd functions (as known, any function defined within the area [−x;x] can be represented as a sum of an even function 0.5∗ [ f (x) +f (−x)] and an odd function 0.5∗ [ f (x) − f (−x)] Korn & Korn
1968) and the odd components of the loading parts of the stress–
strain curves obtained for the three sites are usually larger than theeven components, except the latter intervals in the upper 13 m at PI,
where the stress–strain curves are very sloping (Fig. 3), and the evenand odd components are approximately equal. Even-order higher

harmonics are generated in the latter intervals at PI and at SGK site,where the corresponding stress–strain curves gained noticeable evencomponents, and we can conclude that odd-order non-linearities are
typical for soils, whereas even-order non-linearities appear only incases when the loading (and unloading) parts of stress–strain relations represent functions with noticeable even components.
We should emphasize the difference between the described effects
of higher harmonic generation and the resonant amplification of seismic signals in the upper soft soil layers. As known, resonant spectral
peaks are usually wide and fuzzy, in contrast to spectra shown in
Fig. 4. Here, we consider non-linear distortions of the propagating monochromatic signals, which are caused by the non-linearityof the hysteretic stress–strain relations describing the behaviour ofthe soil layers. For any frequencies of input signals, we will obtainthe same patterns of higher harmonics generation similar to thoseshown in Fig. 4, which do not depend on the boundary conditionsof the models.

C 2005 RAS, GJI, 160, 539–553

Non-linear behaviour during the 1995 Kobe earthquake

545

Figure 4. Velocity spectra of testing monochromatic signals, propagating in soil layers, at depths of locations of the recording devices at PI (0, 16, 32 and83 m), SGK (0, 24.9 and 97 m) and TKS (0, 25 and 100 m), for 10 successive time intervals.


C 2005 RAS, GJI, 160, 539–553

546

O. V. Pavlenko and K. Irikura

The results of testing our models by the Gaussian white noisein 10 successive time intervals are shown in Figs 5(a)–(c) and 6,
and in Tables 1 and 2. Functionals of the Wiener series were calculated by the method of cross-correlation functions (Marmarelis &Marmarelis 1978), and linear, non-linear quadratic, non-linear cubic
components and constant components of the ground response were
estimated according to formulae (2)–(5). Estimates of the zero-,first-, second- and third-order kernels are based on input and output
signals of 500 000 points of duration (more than 40 min). Thesekernels were used to calculate the responses of the linear models
and the non-linear corrections accounting for quadratic and cubic
non-linearities.
Figs 5(a)–(c) and Table 1 show the non-linear components of thesoil responses at PI (Fig. 5a), SGK (Fig. 5b) and TKS (Fig. 5c) sites,
and the non-linear quadratic and cubic components in percentagesof the whole intensity of the responses. At PI, liquefaction occurredin the upper layers and five time intervals are shown in Fig. 5(a),corresponding to different stages of liquefaction. Figs 5(b) and (c)represent examples of the non-linear identification of the soil behaviour at the beginning, middle and final periods of the strong
motion at SGK and TKS sites. The results of the non-linear identification of the soil behaviour in 10 successive time intervals during
the 1995 Kobe earthquake at the three sites are given in Table 1.
A noticeable difference in the shapes of input and output signals
at PI (Fig. 5a) indicates a high non-linearity of the soil behaviourof the upper layers, which increases with developing liquefaction.
A substantial non-linear component appears in the response of soils
at PI even during the first seconds of the strong motion, which increases from ∼40 to ∼60 per cent with developing liquefaction(Fig. 5a, Table 1). At the beginning of the strong motion, the nonlinear cubic component of the response substantially exceeds the
non-linear quadratic one. As liquefaction developed at PI, the nonlinear quadratic component of the response increased, whereas the
non-linear cubic component slightly decreased, so that their partsin the soil response became approximately equal in liquefied soils.
At the same time, the residual part of the soil response, which ispartially the result of higher-order non-linearities, noticeably increased with developing liquefaction, indicating an increase in thehigher-order (the fourth-, fifth-, sixth-order, etc.) non-linearities ofthe soil response at PI. Estimates presented in Table 1 reveal thesetendencies, except the last time interval, where the non-linear components of the soil response are weaker because of the decrease in
the intensity of the strong motion.
The residual parts of the response, remaining after subtractionof linear, non-linear quadratic and non-linear cubic components
from the responses, characterize the parts of the higher-order nonlinearities, as well as inaccuracies of the estimates (Marmarelis &Marmarelis 1978). The inaccuracies are caused by the finite duration and the boundedness of the spectral band of the Gaussian noise,by truncation of tails of the normal distribution of noise amplitudes
and by inaccuracies in calculations. The detailed error analysis andcorrection methods are given in (Marmarelis & Marmarelis 1978).
To evaluate the level of inaccuracies in the obtained estimates, we
performed numerical simulations with linear models, where linear
stress–strain relationships, similar to Hooke’s law, were used insteadof the non-linear hysteretic ones. We found that for the selected
duration of testing Gaussian white noise signals (500 000 points)
estimated non-linear components of the responses of linear models
do not exceed 0.8–1 per cent of the intensity of the response. We
took these values as a rough estimate of the accuracy of our calculations of the Wiener functionals in the non-linear identification.
However, inaccuracies in estimates of the contents of the non-linear

components in the soil response during the 1995 Kobe earthquake
at PI, SGK and TKS sites are evidently higher, up to some per cent,because they account for deviations of the constructed numericalmodels of the soil behaviour during the 1995 Kobe earthquake atthese sites from the real soil behaviour.
At SGK site, changes in the behaviour of the upper layers during
the earthquake were also high and the non-linear part of the response
reached ∼39 per cent, whereas at TKS site it did not exceed ∼13 percent. As seen from Table 1 and Figs 5(b) and (c), the parts of nonlinear components in the soil response at SGK and TKS sites changein accord with the intensity of the strong motion: they are small atthe beginning and at the end of the strong motion, but they growup and reach their maxima in the middle of the strong motion at
moments of its highest intensity. The quadratic, cubic and higherorder non-linear components of the response show similar trends andthe non-linear cubic components are substantially higher than the
non-linear ratic ones. On the whole, the non-linear components of
the response are noticeably higher at SGK than at TKS site (Table 1).
Zero-order kernels h 0 represent quasi-static deformations of thesurface, that is, non-zero values of h 0 indicate slow shifting of thesurface layers to one or the other side in their oscillations. This effectcan be interpreted as a result of accumulation of residual shift deformations on the surface (Zvolinskii 1982). The asymmetry of oscillations is connected with even-order non-linearities (Marmarelis& Marmarelis 1978). Figs 5(a)–(c) show the values of the quasistatic deformations, obtained in our numerical simulations. Higher
values of the quasi-static deformations at PI site and in the second
time interval at SGK site compared with similar values for TKS site
indicate higher contents of the even-order non-linearities in the soil
response at PI and SGK sites. At PI, the quasi-static deformations
increase with developing liquefaction, indicating an increase in theeven-order non-linearities.
Below in Figs 5(a)–(c), values of the first-order Wiener kernelh 1 (τ ), the diagonal values of Wiener kernels h 2 (τ 1 , τ 2 ) and h 3 (τ 1 ,
τ 2 , τ 3 ), and averaged within the upper layers stress–strain relations,which are due to the propagation of the Gaussian white noise in thesoil profiles at PI, SGK and TKS sites, are presented. To discuss
the meaning of the Wiener kernels, let us represent the loading (orunloading) parts of the stress–strain relations as the power series
σ (ε) = a1 ε + a2 ε2 + a3 ε 3 + a4 ε4 + a5 ε5 + . . . ,(6)where σ is stress, ε is strain and a 1 , a 2 , a 3 , a 4 , a 5 are constantcoefficients. The first term in the series a 1 ε describes the linear dependence (the Hooke’s law); the others are non-linear corrections. In
this representation, the first-order kernel h 1 (τ ) is related to the linear
term, so that the maximum of the function h 1 (τ ), corresponding to
the time of seismic wave propagation from the depth of location of
the deepest device to the surface, estimates coefficient a 1 , whereash 2 (τ 1 , τ 2 ) and h 3 (τ 1 , τ 2 , τ 3 ) are related to the quadratic and cubic
non-linear corrections, and maxima or minima of h 2 (τ , τ ) and h 3 (τ ,τ , τ ) [achieved at the same time delays as maxima of h 1 (τ )] evaluate
coefficients a 2 and a 3 , respectively.
We can see from Figs 5(a)–(c) that the first-order kernels h 1 (τ )and the diagonal values of the second- and third-order kernels h 2 (τ 1 ,τ 2 ) and h 3 (τ 1 , τ 2 , τ 3 ) change with time during the strong motion,and the diagonal values of h 3 (τ 1 , τ 2 , τ 3 ) are usually substantiallyhigher than the diagonal values of h 2 (τ 1 , τ 2 ) and h 1 (τ ), especially atthe beginning of the strong motion in the domain of small strains andlarge values of a 1 . The obtained estimates indicate a high variabilityof functions h 2 (τ 1 , τ 2 ) and h 3 (τ 1 , τ 2 , τ 3 ), especially h 3 (τ 1 , τ 2 ,τ 3 ), during the strong motion and a wide range of their variations.

C 2005 RAS, GJI, 160, 539–553

Figure 5. Testing of soil layers with the Gaussian white noise at sites (a) PI (NS component, intervals 1, 3, 5, 7 and 9) for five different stages of liquefaction, (b) SGK (NS component, intervals 1, 2 and 7) and (c)TKS (EW component, intervals 1, 4 and 8). 1 is the input signal, 2 is the output, 3 is the response predicted by a linear model, 4 is the difference between the response of the system and the response predicted by alinear model, 5 is the non-linear correction as a result of quadratic non-linearity predicted by kernel {h 2 }, 6 is the difference between the response of the system and the response predicted by model {h 0 , h 1 , h 2 },7 is the non-linear correction as a result of cubic non-linearity predicted by kernel {h 3 }, and 8 is the difference between the response of the system and the response predicted by model {h 0 , h 1 , h 2 , h 3 }. Below are
shown corresponding average stress–strain relations for groups of layers, represented as even and odd functions, and diagonal value of kernels of the first, second and third orders. Values o/e are the ratios of theamplitudes of the even and odd functions. Zero-order kernels h 0 increase with developing liquefaction.

Non-linear behaviour during the 1995 Kobe earthquakeC 2005 RAS, GJI, 160, 539–553

547

548

O. V. Pavlenko and K. Irikura

Figure 5. (Continued.)

As Figs 5(a)–(c) show, extreme values of h 3 (τ 1 , τ 2 , τ 3 ) can be tens or
hundreds of times higher than h 1 (τ ). Thus, our results explain highvalues of effective parameters of elastic non-linearity in subsurface
soils of ∼102 ÷ 104 , obtained in field experiments (Gushchin &Shalashov 1981; Shalashov 1984; Groshkov et al. 1991), as well as
their large variations. Because seismology deals with a wide range ofstrain levels and, moreover, soils can change their properties during
the strong motion, large variations of the effective parameters of
non-linearity are easy to understand. According to our estimates,
the range of strains during the 1995 Kobe earthquake at the threesites was approximately 10−5 –10−3 .
Loading and unloading parts of the hysteretic stress–strain curves
can be represented as sums of even and odd functions, and our results
confirm that the relationships between the even and odd componentsof the stress–strain curves define the relationships between even-order and odd-order non-linear components of the soil response.
Figs 5(a)–(c) show even and odd components of the loading partsof the stress–strain curves together with the stress–strain curves:
convex-up arcs crossing the hysteretic loops represent the even components and diagonal lines inside the loops show the odd components. Ratios of the amplitudes of the odd to the even components
are designated as o/e and also shown in the figure.
As liquefaction developed at PI, the stress–strain relations in theupper layers became more and more sloping, the even componentsof the loading parts of the stress–strain curves increased and theeven-order non-linear components in the soil response increased(Fig. 5a). Loading parts of the stress–strain relations describing the
behaviour of the upper soil layers at SGK site also gained substantial
even components, when the intensity of the strong motion increased
to its highest values; as a result, noticeable even-order non-linear
C 2005 RAS, GJI, 160, 539–553

Non-linear behaviour during the 1995 Kobe earthquake

549

Figure 5. (Continue.)components appeared in the soil response at SGK (Fig. 5b). At TKS
site, the intensity of the strong motion was not so high, the loading parts of the stress–strain relations describing the soil behaviour
represent predominantly odd functions and even-order non-linear
components in the soil response were weak (Fig. 5c).
At the three sites, temporal changes in the Wiener kernels during
the strong motion are significant (Figs 5a–c). At PI, maxima of
the first-order Wiener kernels h 1 (τ ) decrease and shift to larger
time delays, which is a result of the increased sloping of the stress–
strain relations. Amplification of seismic waves decrease and shear
moduli in the upper layers decrease with developing liquefaction,
whereas the traveltimes of seismic waves increase. At the same time,
the first-order kernels become smoother and wider, which testifies
to the rising predominance of low-frequency oscillations on thesurface, and the second-order kernels h 2 (τ 1 , τ 2 ) grow and smoothC 2005 RAS, GJI, 160, 539–553

with developing liquefaction. At SGK and TKS sites, changes in
the Wiener kernels are also substantial, and we can note that anincrease in the intensity of the strong motion leads to smoothing of
the Wiener kernels and decreasing their amplitudes.
Fig. 6 demonstrates spectra of testing Gaussian white noise signals, obtained in our numerical simulations, at depths of locationsof the recording devices at the three sites in 10 successive timeintervals. These spectra show an important tendency of spectral
transformations in the soil layers, such as, spectra of signals on thesurface tend to take the form of E( f ) ∼ f −k (Fig. 6). Shapes ofspectra at intermediate depths (16 and 32 m at PI; 24.9 m at SGK
and 25 m at TKS) testify to the resonant amplification of some
spectral components in the inner layers. At the same time, interactions of the spectral components take place together with their resonant amplification and lead to generation of combination-frequency

550

O. V. Pavlenko and K. Irikura

Figure 6. Velocity spectra of the testing Gaussian white noise, propagating in soil layers, at depths of locations of the recording devices at PI (0, 16, 32 and83 m), SGK (0, 24.9 and 97 m) and TKS (0, 25 and 100 m), for 10 successive time intervals.


C 2005 RAS, GJI, 160, 539–553

Non-linear behaviour during the 1995 Kobe earthquake

551

Table 1. Non-linear components of the soil response estimated in the numerical simulation with the Gaussian white noise.

Port
Island

EW

NS

SGK

EW

NS

TKS

EW

NS

Number oftime interval

1

2

3

4

5

6

7

8

9

10

nl
nl-2
nl-3
res
nl
nl-2
nl-3
res

41.9
8.7
14.3
21.3
39.1
7.2
12.9
20.3

49.7
9.4
13.2
27.3
51.2
9.5
14
27.5

37.7
7.3
11.5
20.2
40.2
7.3
13.6
20.8

54.5
10.6
12.9
30.5
39.3
7.3
13.4
20.3

46
8.7
10.6
28.6
49.4
9.2
12.4
28.1

53.1
10.1
9.2
31.3
44.8
8.3
13.2
23.8

53
10.2
9.0
30.6
59
10.5
9.9
34.9

54.2
10.5
9.1
31.3
53.2
9.3
9.5
31.5

59.3
11.8
9.9
33.3
56.9
10.7
9.1
33.4

48.3
8.6
8.7
28.8
47.2
9.1
9
27.3

nl
nl-2
nl-3
res
nl
nl-2
nl-3
res

4.3
0.8
2.3
2.2
3.4
0.6
1.7
1.8

39.1
7.3
12.9
20.6
27.4
5.1
9.2
13.8

35.2
6.7
11.9
17.7
37.5
6.6
18.8
26.4

28.8
5.3
11.2
14.9
33.9
6.4
12.2
17.3

33.6
6.4
11.4
17.7
14.8
2.8
6.8
7.4

11.1
2.1
4.9
5.4
22.5
4.4
8.4
11.2

10.3
2
4.9
4.9
11.3
2.1
5
5.5

9.7
1.8
4.4
4.7
8.2
1.6
3.7
4.1

8.9
1.7
3.8
4.6
7.6
1.4
3.6
4

8.9
1.7
3.7
4.5
6.2
1.2
2.5
3.1

nl
nl-2
nl-3
res
nl
nl-2
nl-3
res

1.8
0.3
0.8
1.2
1.4
0.3
0.6
1.1

4.3
0.8
2.2
2.3
3.5
0.7
2
2.2

6.2
1.2
3
2.8
9.1
1.7
5.2
4.4

7.3
1.4
3.7
3.3
12.7
2.4
6.6
5.7

6.1
1.2
2.9
2.8
10.5
2
5.3
4.6

7.5
1.4
4.4
3.7
11.1
2.1
5.8
5.1

5.7
1.1
3
2.8
9.1
1.7
5.1
4.5

4.6
0.9
2.3
2.4
3.8
0.7
2.3
2.4

2.9
0.5
1.3
1.6
4.1
0.8
2
2

3.2
0.6
1.9
2
2.5
0.5
1.1
1.5

nl is the non-linear part of the response, which represents the deviation of the response of a real system from the response of the linear model constructed by
the zero- and first-order kernels {h 0 , h 1 };nl-2 is the non-linear correction as a result of quadratic non-linearity predicted by kernel {h 2 };nl-3 is the non-linear correction as a result of cubic non-linearity predicted by kernel {h 3 };res- the deviation of the system response from the response of the non-linear model predicted by kernels {h 0 , h 1 , h 2 , h 3 }.
Table 2. Coefficient k estimated in the numerical simulation with the Gaussian white noise in 10 successive time intervals and from the surface velosigrammsof the 1995 Kobe earthquake (the last column).
Number oftime interval

1

2

3

4

5

6

7

8

9

10

EQ

Port
Island

EW
NS

−1.33
−1.31

−1.3
−1.31

−1.35
−1.35

−1.18
−1.37

−1.24
−1.28

−1.23
−1.43

−1.31
−1.25

−1.13
−1.26

−1.28
−1.22

−1.32
−1.44

−1.21
−1.39

SGR

EW
NS

−0.31
−0.25

−1.28
−1.41

−1.51
−1.26

−1.1
−1.51

−1.43
−1.48

−1.46
−1.48

−1.3
−1.45

−1.28
−1.43

−1.08
−1.42

−1.08
−1.34

−1.46
−1.39

TKS

EW
NS

−0.32
−0.29

−0.51
−0.45

−0.62
−0.79

−0.69
−0.95

−0.63
−0.86

−0.69
−0.89

−0.6
−0.79

−0.52
−0.47

−0.4
−0.49

−0.43
−0.38

−1.45
−1.56

harmonics in the low and high frequency ranges, amplitudes ofcombination harmonics being related to their frequencies (Zarembo& Krasil’nikov 1966). As a result, the energy is redistributed overthe spectral band so that low-frequency components increase, sharp
spectral peaks disappear and the resulting spectrum of the output
tends to take the limiting smooth form of E( f ) ∼ f −k (Fig. 6).
Table 2 represents estimates of coefficient k obtained in our numerical simulation with the Gaussian white noise in 10 successive
time intervals and coefficient k estimated from the surface velosigramms of the 1995 Kobe earthquake. As is seen from the table,
at PI, coefficient k does not change much with developing liquefaction and varies within the limits of ∼1.13 to ∼1.44, whereas at
SGK and TKS, variations of coefficient k are substantial. Note thatin cases when spectra of signals on the surface can be accurately
approximated by expression E ∼ f −k , k values are close to 1.3 ÷
C 2005 RAS, GJI, 160, 539–553

1.4 (Fig. 6, Table 2). These cases are marked in bold in Table 2.
Non-linear distortions of the propagating signals were high at PI
and all the spectra of signals on the surface took the smooth form ofE ∼ f −k (Fig. 6). At SGK and TKS sites, non-linearity was not sostrong and spectra of signals on the surface approach their limiting
form only at SGK in the intervals of the highest intensity of thestrong motion, and in these cases coefficients k take values similar
to those at PI, i.e. 1.26–1.41. At TKS site, spectra of signals on thesurface can not be accurately approximated by expression E ∼ f −k
and coefficients k vary during the strong motion, increasing in its
middle period.
Values of coefficient k estimated from the surface velosigrammsof the 1995 Kobe earthquake are shown in the last column in Table 2.
These values are close to those obtained in the numerical simulation
with the Gaussian white noise, especially for PI. Velocity spectra of

552

O. V. Pavlenko and K. Irikura

Figure 7. Acceleration spectra of the main shock of the 1995 Kobe earthquake at PI, SGK and TKS sites.

the 1995 Kobe earthquake are more close to the form E ∼ f −k inthe case of PI and less close in the case of TKS site.
Thus, intense seismic waves can suffer substantial distortions
in subsurface soils; as a result, the information about the spectral
composition of input signals is lost. However, the limiting shapeof spectrum can be attained only in numerical simulation, when
we analyse long-duration signals. In practice, we observe only thetendency of increasing low-frequency components and smoothing
spectral peaks, because spectra are calculated for rather short intervals of records. Acceleration spectra of the 1995 Kobe earthquake
are shown in Fig. 7 and the tendency is clearly seen: upward propagating seismic signals gain in the low-frequency components andloose sharp spectral peaks at high frequencies. As a whole, acceleration spectra of signals on the surface and at depths possess few
common features.

4 C O N C LU S I O N S
We carried out the non-linear identification and estimated linear andnon-linear components of the soil response during the 1995 Kobe
earthquake at PI, SGK and TKS sites. The conclusions of the study
are as follows.
The non-linear components of the soil response decreased withthe distance from the fault plane from values of ∼60 per cent at PI (2km from the fault plane) to ∼ 40 per cent at SGK (6 km from the fault
plane) and to ∼13 per cent at TKS (15 km from the fault plane; thestrain range is 10−5 –10−3 ). Odd-order non-linear components usually predominate in the soil response, whereas even-order non-linear
components become comparable to odd-order ones only in cases ofliquefied soils or high intensity of the strong motion, when the loading parts of stress–strain relations in the upper layers gain noticeable
even components. Constant components of seismic wavefields, i.e.
quasi-static deformations of the surface, change in accord with theeven-order non-linear components of the soil response. Liquefaction induces a substantial increase in the non-linear components ofthe soil response, from ∼40 to ∼60 per cent of the intensity of the
response at PI, and changes in the predominant types of the soil nonlinearity, which are the result of the rise of the even-order non-linear
components of the response.
The ratio of the even-order non-linear components of the soil
response during the strong motion to the odd-order non-linear components corresponds to the ratio of the even to the odd components of
the loading parts of the non-linear hysteretic stress–strain relationsof the upper most non-linear soil layers.
In cases of strong non-linearity, i.e. thick sedimentary layers andstrong ground motion, spectra of seismic signals on the surface
tend to take the form of E ∼ f −k . This limiting spectral shape wasachieved at PI and SGK sites during the 1995 Kobe earthquake andk ≈ 1.2 ÷ 1.4 for velocity spectra on the surface.
The proposed methods for processing vertical array records allow extracting maximum information on the soil behaviour from therecords, and help in understanding the mechanisms of transformations of seismic waves in subsurface soil layers and predicting thesoil behaviour during future earthquakes.
AC K N OW L E D G M E N T S
We thank Yoshio Soeda, Hideki Tamai and Naoyuki Nakatsu (Kansai
Electric Power Company) for acceleration records of the verticalarrays at SGK and TKS recording sites. The acceleration records ofthe vertical array in PI were kindly provided by Kobe City. This study
was partially supported by Grant-in-Aid for Scientific Research, Nos
11792026 and 00099 Japanese Society for the promotion of Science,from the Ministry of Education, Science, Sports and Culture, Japan.

REFERENCES
Aleshin, A.S., Gushchin, V.V., Krekov, M.M., Nikolaev, A.V. & Shalashov,G.M., 1981. Experimental studies of nonlinear interaction of seismic surface waves, Dokl. Akad. Nauk SSSR, 260(3), 574–575.

C 2005 RAS, GJI, 160, 539–553

Non-linear behaviour during the 1995 Kobe earthquake
Archuleta, R.J., 1998. Direct observations of nonlinearity in accelerogramms, in The Effects of Surface Geology on Seismic Motion, Vol. 1,
pp. 787–792, eds Irikura, K., Okada, H. and Sasatani, T., Balkema, Rotterdam.
Archuleta, R.J. & Steidl, J.H., 2001. Borehole array data: open windows into
site response. In: 10th International Conference on Soil Dynamics andEarthquake Engineering, Philadelphia, USA, October 7–10, 2001, Volume of Extended Abstracts, vi, ed. Zerba, A., Drexel University, Philadelphia.
Bard, P.-Y. & Pitilakis, K., 1995. Seismic zonation and ground motion interface. In: Proc. of the 5th International conference on Seismic Zonation,
October. 17–19, Nice, France, Vol. III, pp. 2127–2152, eds Pecker, A.,
Cluff, L., Ouest Editions, Versailles, France.
Dimitriu, P.P., 1988. Self-modulation and recurrence phenomena in vibratorinduced steady-state sinusoidal ground vibrations, Phys. Earth planet. Int.,50, 74–82.
Dimitriu, P.P., 1990. Preliminary results of vibrator-aided experiments in
non-linear seismology conducted at Uetze, F.R.G., Phys. Earth planet.
Int., 63, 172–180.
Field, E.H., Johnson, P.A., Beresnev, I.A. & Zeng, Y., 1997. Nonlinear
ground-motion amplification by sediments during the 1994 Northridge
earthquake, Nature, 390, 599–602.
Groshkov, A.L. & Shalashov, G.M., 1986. Equations of nonlinear dynamics of elastic-relaxation media, Dokl. Akad. Nauk SSSR, 290(4), 825–
827.
Groshkov, A.L., Kalimulin, R.R. & Shalashov, G.M., 1991. Cubic nonlinear
seismic effects, Dokl. Akad. Nauk SSSR, 315(1), 65–68.
Gushchin, V.V. & Shalashov, G.M., 1981. On the possibility of applicationof nonlinear seismic effects to problems of vibro-radiation of the earth, inStudy of the Earth by nonexplosion seismic sources, Nauka, Moscow, pp.
144–155.
Hardin, B.O. & Drnevich, V.P., 1972. Shear modulus and damping in soils:
Design equations and curves, J. Soil Mech. Found. Div., 98, 667–692.
Joyner, W.B. & Chen, T.F., 1975. Calculation of nonlinear ground responsein earthquakes, Bull. seism. Soc. Am., 65, 1315–1336.
Kadomtsev, B.B. & Karpman, V.I., 1871. Nonlinear waves, Uspekhi fiz. nauk,
103(2), 27–48.
Kawase, H., Satoh, T., Fukutake, K. & Irikura, K., 1995. Borehole records
observed at the Port Island in Kobe during the Hyogo-ken Nanbu earthquake of 1995 and its simulation, J. Constr. Eng., AIJ, 475, 83–92.
Kazama, M., Yamaguchi, A. & Yanagisawa, E., 1998. Seismic behaviour ofC 2005 RAS, GJI, 160, 539–553

553

an underlying alluvial clay on man-made islands during the HyogokenNanbu earthquake, Soils and Foundations, Special Issue, 2, 23–32.
Korn, G.A. & Korn, T.M., 1968. Mathematical handbook for scientists
and engineers, McGraw-Hill Book Company New York, San Francisco,
Toronto, London, Sydney.
Kudo K., 1995. Practical Estimates of Site Response State-of-the-Art Report.
In: Proc. of the 5th International conference on Seismic Zonation, October
17–19, Nice, France, Ouest Editions, Versailles, France, Vol. III, 1103–
1112.
Lund, F., 1983. Interpretation of the precursor to the 1960 Great Chilean
Earthquake as a seismic solitary wave, Pageoph, 121(1), 17–26.
Marmarelis, P.Z. & Marmarelis, V.Z., 1978. Analysis of Physiological Systems, The White-Noise Approach, Plenum Press, New York and London.
O’Connell, D.R.H., 1999. Replication of apparent nonlinear seismic response with linear wave propagation models, Science, 283, 2045–2050.
Pavlenko, O.V., 2001. Nonlinear Seismic Effects in Soils: Numerical Simulation and Study, Bull. seism. Soc. Am., 91, 381–396.
Pavlenko, O.V. & Irikura, K., 2002. Types of elastic nonlinearity of sedimentary soils, Geophys. Res. Lett., 29(19), 36–36-4.
Pavlenko, O.V. & Irikura, K., 2003. Estimation of nonlinear time-dependent
soil behaviour in strong ground motion based on vertical array data, Pureappl. Geophys., accepted.
Sato, K., Kokusho, T., Matsumoto, M. & Yamada, E., 1996. Nonlinear seismic response and soil property during 1995 Hyogoken Nanbu earthquakes
Soils and Foundations, Special issue on Geotechnical Aspects of the January 17, 1995 Hyogoken Nanbu earthquake, 41–52.
Shalashov, G.M., 1984. Cross-modulation of acoustic waves on the cubic
nonlinearity in solids, Acoustical Journal, XXX(3), 386–390.
Van Den Abeele, K.E.-A. & Johnson, P.A., 1996. Elastic pulsed wave propagation in media with second- or higher-order nonlinearity. Part II. Simulation of experimental measurements on Berea sandstone, J. acoust. Soc.
Am., 99, 3346–3352.
Vasil’ev, Yu.I., Ivanova, L.A. & Shcherbo, M.N., 1969. Measurement ofstress and strain in a soil during the propagation of waves from explosions,Izv. Akad. Nauk SSSR, Fiz. Zemli, 1, 21–37.
Zarembo, L.K. & Krasil’nikov, V.A., 1966. Introduction to Nonlinear Acoustics, Nauka, Moscow.
Zimenkov, S.V. & Nazarov, V.E., 1993. Nonlinear acoustic effects in rocksamples, Izv. Akad. Nauk SSSR, Fiz. Zemli, 1, 13–18.
Zvolinskii, N.V., 1982. Wave processes in non-elastic media, Problems ofengineering seismology, 23, 4–19.