Ó Birkhäuser Verlag, Basel, 2003

Pure appl. geophys. 160 (2003) 2365–2379
0033 – 4553/03/122365 – 15
DOI 10.1007/s00024-003-2398-9

Pure and Applied Geophysics

Estimation of Nonlinear Time-dependent Soil Behavior in
Strong Ground Motion Based on Vertical Array Data
O. V. PAVLENKO1 and K. IRIKURA2

Abstract — To improve our understanding of nonlinear elastic properties of soils, a method is
proposed of estimation of stress-strain relations of soils in situ in strong ground motion based on verticalarray data. Strong motion records provided by seismic vertical arrays allow estimation of nonlinear stressstrain relations in soil layers at diﬀerent depths, from the surface down to the location of the deepest
device. As an example, records obtained during the main shock of the 1995 Kobe earthquake at PortIsland, SGK, and TKS sites were used to estimate the stress-strain relations in the soil proﬁles. For
diﬀerent layers, diﬀerent types of nonlinear stress-strain relations were selected, according to the proﬁling
data. To account for temporal changes in the soil behavior, consecutive parts of records were examined,
and for successive time intervals, the relations were found showing the best-ﬁt approximation to theobserved data. At Port Island and SGK sites, where the strongest accelerations were recorded, the obtained
stress-strain relations showed systematic changes in the upper layers (0–14 m), such as, a progressive
reduction of the slopes of the stress-strain curves due to liquefaction at Port Island and reduction andrecovery of the slopes at SGK and TKS sites. At the three sites, the stress-strain relations remained stablein layers below 11–14 m. Thus, the proposed approach gives us a representation of the soil behavior inlayers at diﬀerent depths in strong ground motion; it allows calculation of the propagation of arbitrary
seismic signals in the studied proﬁles and estimation of nonlinear components in the ground response by
the nonlinear system identiﬁcation technique. The method can also be applied to evaluate the groundresponse at sites where proﬁling data are available and an imposed motion can be estimated.
Key words: Strong ground motion, nonlinear soil behavior, liquefaction, stress-strain relations,seismic vertical arrays.

Introduction
Recent earthquakes, such as the 1985 Michoacan earthquake, the 1989 Loma
Prieta earthquake, the 1994 Northridge earthquake, and the 1995 Kobe earthquake
provided new experimental data on the soil behavior in strong ground motion, inparticular, on the liquefaction phenomena, and discussions on the nonlinearity of soil
behavior were induced (LOMNITZ et al., 1995; AGUIRRE and IRIKURA, 1997; FIELD
et al., 1997; O’CONNELL, 1999, etc.). Though nonlinear elastic properties of soils were1
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

2366

O. V. Pavlenko and K. Irikura

Pure appl. geophys.,studied in multiple laboratory experiments, and valuable laboratory experimental
data are accumulated sometimes, this is not suﬃcient for understanding the soil
behavior in situ, because soils often represent multicomponent systems containing
water, air, gases, etc., and strong ground motion can induce movement andredistribution of these components, i.e., changes in the properties of the soils.
Experimental data on the soil behavior in strong motion in situ are still few,
fragmental, and non-representative; and accumulation of these data is important for
improving our understanding of soil behavior in strong motion.
In strong ground motion Hooke’s law does not hold for subsurface soils, i.e., soils
should be taken as nonlinear systems transforming incident seismic signals into
movement on the surface. For studying nonlinear properties of systems, eﬀective
methods are developed in system analysis, so-called nonlinear system identiﬁcation
technique, based on the determination of higher-order impulse characteristics of thesystems. An output of a nonlinear system is represented as the Volterra-Wiener series,
i.e., a sum of the response of a linear system to the input signal and a number ofnonlinear corrections, which are due to quadratic, cubic nonlinearity, and nonlinearities of higher (the 4-th, 5-th, etc.) orders. If we know the input and output of anonlinear system, we can judge regarding the types and quantitative characteristics ofthe system nonlinearity (MARMARELIS and MARMARELIS, 1978). Nonlinear identiﬁcation of soils in various geotechnical conditions seems to be promising, because it
allows a better understanding of the behavior of soils and structures in strong groundmotion. However, to apply methods of system analysis to studying nonlinear
properties of soils, knowledge of stress-strain relations in the soil layers in strong
motions is required. In this paper, a method of estimation of nonlinear stress-strain
relations in soils in strong ground motion is proposed based on vertical array data.
Numerous methods and programs developed for calculating the groundresponse in strong motion in various conditions do not allow estimation ofstress-strain relations in soil layers in situ. Moreover, in cases of strong
nonlinearity, there often remains some disagreement between the observed andsimulated records. As is known, equivalent linear models (SHAKE, QUAD-4,
FEADAM, LUSH, FLUSH) are not applicable for calculation of such complex
phenomena as soil liquefaction. Programs DESRA (LEE and FINN, 1978), TARA
and their modiﬁcations (FINN et al., 1986; FINN and YOGENDRAKUMAR, 1989)
allow determination of the possible level of the pore pressure and the possibility ofliquefaction, and they can be applied for the analysis of soil behavior afterliquefaction. Changes in the pore pressure are related to the volumetric deformations in soils in drained conditions, and one-dimensional diﬀusion is included in thealgorithm. Programs DYSAC2, DYNAFLOW, and SWANDYNE are considered
to provide the best results (ARULANANDAN et al., 1995). Equations of motion of
the liquid and solid phases are related to the equation of conservation of matter.
Generation and dissipation of the pore pressure are connected with the deformation of the solid matrix due to the Biot equations (BIOT, 1956). However, in every

Vol. 160, 2003

Estimation of Nonlinear Soil Behavior

2367

case simpliﬁcations and assumptions are applied, concerning the medium properties, as well as the mechanisms of the processes, therefore, any uncertainties andmistakes in modelling lead to an improper calculation of the soil movement.
At the same time, records of strong ground motion provided by seismic verticalarrays allow estimation of stress-strain relations in soil layers in situ. This paper
describes the method of estimation of stress-strain relations and its application to the1995 Kobe earthquake. The method allows us to trace temporal changes in the stressstrain relations. Since the estimates are based only on real measurements, they are freeof theoretical approximations and physical assumptions concerning mechanisms ofprocesses arising in the medium in strong ground motion.

Data and Method
Data were processed for three recording sites, Port-Island, SGK, and TKS.
Distances to the closest point on the fault line are 2 km, 6 km, and 16 km, respectively.
Figure 1, taken from the paper by SATO et al. (1996), shows the locations of the sites,
the major principal axes, and the epicentres of the main shock and aftershocks
summarised by Disaster Prevention Research Institute of Kyoto University.
At Port Island, the vertical array contains four three-component accelerometers,
installed at GL-0 m, GL-16 m, GL-32 m, and GL-83 m; the arrays at SGK and
TKS sites consist of three three-component devices at GL-0 m, GL-24.9 m, andGL-97 m, and GL-0 m, GL-25 m, and GL-100 m, respectively (Fig. 2). We checked
the directional drifts of the accelerometers by calculating the horizontal orbits of
the long-period particle motions for the main shock and the aftershocks at diﬀerent
depths at the three sites. At Port Island, N19°W rotation at GL-83 m was detected
and corrected; at SGK site, a reverse of NS component and N6°W rotation at GL24.9 m and a reverse of NS component and N34°E rotation at GL-97 m weredetected and corrected. At TKS site, N23°W rotation at GL-25 m and N9°E
rotation at GL-100 m were detected and corrected. All these corrections agree withthe conclusions of other authors.
The material at the three sites are similar to one another: reclaimed soil, clays, sands,
and gravel (Fig. 2). The proﬁling data at Port Island used for nonlinear simulation weretaken from (AGUIRRE and IRIKURA, 1997). Shear wave velocity, shear modulus
degradation, and maximum shear strain at SGK and TKS sites were taken from SOEDA
et al. (1999), These data were used to calculate the shear stress in failure smax and thelow-strain shear modulus Gmax in diﬀerent layers, following SEED et al. (1984) and SUN
et al. (1988). P-wave velocities, when exceeding 1000 m/s, indicate saturation of soils
with water. The underground water level lay at about 13 m at Port Island and 3–5 m
at SGK and TKS. The deepest layers were Pleistocene gravely soils and the upper layers
consist of alternating Pleistocene gravel/clay layers, Holocene sand/clay layers and ﬁll(SATO et al., 1996).

2368

O. V. Pavlenko and K. Irikura

Pure appl. geophys.,

Figure 1
Locations of vertical array sites Port Island (PI), SGK, and TKS around Osaka bay, the major principal
axes, and the epicentres of the main shock and aftershocks (derived from SATO et al., 1996).

To estimate the stress-strain relations in diﬀerent layers, vertical array records ofstrong motion and the proﬁling data were used. We applied the modiﬁed algorithmof calculation of the propagation of a vertically incident shear wave in soil layers up
to the surface, based on the lumped mass method (JOYNER and CHEN, 1975), andcombined it with the method of ‘‘trial and error’’ to ﬁnd the stress-strain relations
showing the best agreement between the observed and simulated records. In our
computations, stress and strain are normalised in the manner used by HARDIN and
DRNEVICH (1972): stress is normalised by multiplying by 1/smax, and strain isnormalised by multiplying by Gmax/smax. For calculations, the studied medium fromthe surface down to the location of the deepest device was divided into groups oflayers, for which certain types of stress-strain relations were assumed (Fig. 2). Threemain types of stress-strain curves were considered:(1) Those that are similar to laboratory experiments by HARDIN and DRNEVICH(1972), to describe the behavior of dense soils at depths;(2) Those of ‘‘soft’’ type, similar to type (1), but with greater slope, being close tohorizontal for large strains, for liqueﬁed soils;

Figure 2
The ground proﬁles and characteristic types of stress-strain relations at diﬀerent depths at Port Island, SGK, and TKS sites.

Vol. 160, 2003
Estimation of Nonlinear Soil Behavior
2369

2370

O. V. Pavlenko and K. Irikura

Pure appl. geophys.,(3) Those of ‘‘hard’’ type, declining to the stress axis at large strains, for watersaturated soils (terms ‘‘soft’’ and ‘‘hard’’ type stress-strain curves were introduced by
ZVOLINSKII, 1982).
Sets of such curves were generated, and item-by-item examination was applied toﬁnd groups of curves showing the best-ﬁt approximation to the observed data.
To account for temporal changes in the soil behavior, the records were divided
into intervals of 1.5-seconds duration. Within each interval, the stress-strain
relations were assumed to be stationary, and vary for diﬀerent intervals.
Calculations were performed (i.e., the ‘‘best-ﬁt’’ stress-strain relations in the layers
were determined) successively, interval by interval. No discontinuities occur at theboundaries of the intervals, because, in the next interval, the whole cycle of loading(or unloading) was recalculated for the new curve from its beginning. Only
horizontal components of records were processed.
Results and Discussion
The results of the simulation with the ‘‘best-ﬁt’’ stress-strain relations show agood agreement between the observed and simulated records at the three sites.
The most intense and complex movements and liquefaction were observed at Port
Island. The strongest horizontal accelerations were measured at the deepest point of 83m, such as, 500 Gal for the EW component and 600 Gal for the NS component.
According to the material structure and the results obtained in the previous analysis(KAWASE et al., 1995), the medium was divided into ﬁve groups of layers, for which
certain types of stress-strain relations were assumed (Fig. 2). Since liquefaction
occurred in the upper layers, curves of type (2) were applied for surface soils above thelevel of the underground water (0–13 m). Soil layers at depths 13–18 m and 27–32.5 m
are reclaimed gravely ﬁll and alluvial sand saturated with water, therefore, curves oftype (3) were used to describe the soil behavior in these layers (Fig. 2). Peaks, or sharp
increases in the amplitudes of the observed records at depths of 16 m and 32 m conﬁrm
that curves of type (3) are relevant in this case. Stress-strain curves of type (1), which areconsidered to be typical for dense soils at depth, were assumed for alluvial clay layers at
18–27 m and for deep diluvial layers below 32.5 m.
For successive 1.5-second time intervals, the groups of stress-strain relations were
found, showing the best-ﬁt approximation to the observed records. Figure 3 shows
the observed and simulated records (a) and the obtained stress-strain relations at
diﬀerent depths (b) for ten successive time intervals.
cFigure 3
The acceleration time history of the main shock in Port Island, observed and simulated (a), and theobtained stress-strain relations changing with time (b). The axes scales of the stress-strain relations are in
relative units; the same for all time intervals at a given depth.

Vol. 160, 2003

Estimation of Nonlinear Soil Behavior

2371

2372

O. V. Pavlenko and K. Irikura

Figure 3(Contd.)Pure appl. geophys.,Vol. 160, 2003

Estimation of Nonlinear Soil Behavior

2373

The most noticeable changes in the soil behavior were observed in layers near thesurface, in the upper 13 m: the stress-strain relations become more and more slopingwith time, showing a substantial progressive reduction of the shear modulus andliquefaction. At depths 13–18 m and 27–32.5 m, the obtained stress-strain curves
show a slight reduction and the following recovery of the shear modulus. Below
32.5 m, no changes in the soil behavior are observed (Fig. 3b).
At SGK and TKS sites no liquefaction occurred, though sand boils observed
after the quake around the TKS site indicate that liquefaction took place in thevicinities of this site.
Maximum accelerations recorded at SGK site were also high (Fig. 4a), up to 650
Gal on the surface for the EW-component. At this site, the soil behavior can be
described by a relatively simple model: the soil proﬁle is divided into two groups oflayers, such as the near-surface alluvial layers (0–11 m), and the deeper layers (below
11 m), which are mostly gravel and clays. High values of P-wave velocities in theupper layers indicate the presence of water (Fig. 2), therefore, curves of type (3) areappropriate for the upper layers. Curves of type (1) were used for the deeper layers.
Figure 4 represents the observed and simulated accelerograms (a) and the obtained
stress-strain relations (b) for SGK site. The behavior of the layers below 11 m isstationary, whereas the stress-strain relations describing the behavior of the upper
layers change with time: slopes of the curves decrease, then increase again, indicating
reduction and recovery of the shear modulus in the upper layers.
At TKS site, maximum recorded accelerations were about 200 Gal, and the soil
proﬁle is represented by a water-saturated reclaimed ﬁll and an alluvial layer in theupper part (0–14 m), and by gravel, clays, and sands in the deeper parts. Since the soil
conditions are similar to that at SGK site, the choice of curves was also similar:
curves of type (3) were selected for the upper 14 m, and curves of type (1) were chosen
for the deeper parts. The observed and simulated accelerograms and the obtained
stress-strain relations are shown in Figures 5a,b. The most intense movements took
place in the layers below the level of the underground water, at depths 4–7 m. At
TKS site, the soil behavior in the upper layers changes with time similar to SGK site:
we observe reduction and the following recovery of the shear modulus, whereas, thebehavior of layers below 14 m is stationary.
The obtained stress-strain relations seem to represent a fairly good approximation
to reality. On one hand, they show a good agreement between the observed andsimulated data. On the other hand, they give a description of the process which is
physically correct, i.e., the stress-strain relations obtained for Port Island show
progressive liquefaction in the upper layers and a stable behavior in the deeper parts.
Reduction and a following recovery of the shear modulus is obtained for the upper
layers at SGK and TKS sites. Similarity of the stress-strain relations obtained for two
horizontal components is an additional factor testifying to the validity of the solution.
As shown above, the choice of the types of stress-strain relations in layers was
physically justiﬁed. To check the possibility of diﬀerent representations (i.e., vertical

O. V. Pavlenko and K. Irikura

Figure 4
The acceleration time history of the main shock at SGK site, observed and simulated (a), and the obtained stress-strain relations changing with time (b). The
axes scales of the stress-strain relations are in relative units; the same for all time intervals at a given depth.

2374
Pure appl. geophys.,Vol. 160, 2003

Estimation of Nonlinear Soil Behavior

2375

Figure 4(Contd.)distributions of the stress-strain relations), calculations were performed several times
for Port Island and SGK sites, with variations in the shapes of curves, however, every
time similar distributions of the stress-strain relations in layers were obtained.
Characteristic features of the observed accelerograms, such as, sharp peaks at GL16 m and GL-32 m at Port Island and GL-0 m at SGK indicated deﬁnitely certain
stress-strain relations for the corresponding time intervals. At Port Island, the soil
behavior was the most complicated, and the obtained vertical distribution of thestress-strain relations agrees well with the results of KAWASE et al. (1995), wherestress-strain relations in diﬀerent layers were determined by numerical simulation
based on the pore water pressure measurements and the Biot equations. This also cantestify to the validity of our solution.
All things considered, for the given proﬁling data, the obtained distributions ofstress-strain relations in layers at Port Island, SGK, and TKS sites seem to be good
approximations to reality. The number of the recording devices is an importantparameter, determining the accuracy of the description of the system. The four threecomponent accelerometers at Port Island provide better information than the threeO. V. Pavlenko and K. Irikura

Figure 5
The acceleration time history of the main shock at TKS site, observed and simulated (a), and the obtained stress-strain relations (b). The axes scales of thestress-strain relations are in relative units; the same for all time intervals at a given depth.

2376
Pure appl. geophys.,Vol. 160, 2003

Estimation of Nonlinear Soil Behavior

2377

Figure 5(Contd.)devices at SGK and TKS, and allow a more detailed description of the soil behavior
in the layers.
Conclusions
Thus, time-dependent models of nonlinear soil behavior in the 1995 Kobe
earthquake are constructed for three vertical array recording sites at Port Island,
SGK, and TKS. The obtained stress-strain relations indicate temporal changes in the2378

O. V. Pavlenko and K. Irikura

Pure appl. geophys.,soil behavior in the upper layers and a stationarity in the deeper parts. At Port
Island, where the strongest horizontal accelerations were measured (700 Gal at
83 m and 340 Gal on the surface), the most complicated vertical distribution ofstress-strain relations was observed.
For weaker motions at SGK (300 Gal at 97 m and 650 Gal on the surface)
and TKS (100 Gal at 100 m and 200 Gal on the surface), simpler behavior wasobserved. The medium structures are similar at these three sites, therefore, diﬀerencesin the soil behavior are due to the diﬀerences in the acceleration amplitudes of the
input signals. The soil behavior shows increased complexity as amplitudes become
greater. Underground water plays an important role, since it inﬂuences soil behavior
in strong motion and the possibility of liquefaction.
Our results conﬁrm that in strong ground motion, soils not only change
the dynamic parameters of the propagating seismic waves, but they also
substantially change their properties, and these phenomena are inseparable fromeach other.
Our proposed method of estimating the stress-strain relations in soil layers atvertical array sites improves our understanding of processes arising in soils subjected
to strong motion. Knowledge of stress-strain relations in soil layers allows us toperform further studies by calculating the propagation of seismic signals of various
types in the soil proﬁles. We can calculate the propagation of the varieties of imposed
motions—to estimate a range of possible ground responses; the propagation of
monochromatic signals—to investigate the generation of higher harmonics, and the
propagation of the Gaussian white noise—to apply the nonlinear system identiﬁcation technique and determine the types and quantitative characteristics of thenonlinearity of the system of soil layers.
Since the selected ‘‘best-ﬁt’’ curves are determined by the proﬁling data, themethod allows constructing models of soil behavior in strong motions invarious geotechnical conditions and gives us knowledge, which can be usefulin predicting ground motion in future earthquakes. The method can be appliedin any site where proﬁling data are available and an imposed motion can be
estimated.

Acknowledgements
We thank Yoshio Soeda, Hideki Tamai, and Naoyuki Nakatsu (Kansai ElectricPower Company) for acceleration records of the vertical arrays at SGK and TKS
recording sites. The acceleration records of the vertical array in Port Island werekindly provided by Kobe City. This study was partially supported by Grant-in-Aid
for Scientiﬁc Research, No.11792026 and No.00099 (JSPS), from the Ministry of
Education, Science, Sports, and Culture, Japan.

Vol. 160, 2003

Estimation of Nonlinear Soil Behavior

2379

REFERENCESAGUIRRE, J. and IRIKURA, K. (1997), Nonlinearity, Liquefaction, and Velocity Variation of Soft Soil Layers
in Port Island, Kobe, during the Hyogo-ken Nanbu Earthquake, Bull. Seismol. Soc. Am. 87, 5, 1244–1258.
ARULANANDAN, K., MANZARI, M., ZENG, X., FAGAN, M., SCOTT, R. F. and TAN, T. S., Signiﬁcance of the
VELACS project to the solution of boundary value problems in geotechnical engineering. In Third Inter.
Conf. on Recent Advances in Geotechnical Earthq. Engin. and Soil Dyn. (St. Louis, MO, 2 1995), 825–832.
BIOT, M. A. (1956), Theory of Propagation of Elastic Waves in a Fluid-saturated Porous Solid, J. Acoust.
Soc. Am. 28, 2, 168–191.
FIELD, E. H., JOHNSON, P. A., BERESNEV, I. A. and ZENG, Y. (1997), Nonlinear Ground-motion
Ampliﬁcation by Sediments during the 1994 Northridge Earthquake, Nature 390, 6660, 599–602.
FINN, W. D. L., YOGENDRAKUMAR, M., YOSHIDA, N. and YOSHIDA, M., TARA-3: A Program for
Nonlinear Static and Dynamic Eﬀective Stress Analysis (Soil Dynamics Group, Univ. of British
Columbia, Vancouver 1986).
FINN, W. D. L. and YOGENDRAKUMAR, M., TARA-3FL—program for Analysis of Liquefaction Induced
Flow Deformations (Department of Civil Engineering, Univ. of British Columbia, Vancouver 1989).
HARDIN, B. O. and DRNEVICH, V. P. (1972), Shear Modulus and Damping in Soils: Design Equations andCurves, Proc. Am. Soc. Civil Eng., J. Soil Mech. Found. Div. 98, 667–692.
JOYNER, W. B. and CHEN, T. F. (1975), Calculation of Nonlinear Ground Response in Earthquakes, Bull.
Seismol. Soc. Am. 65, 1315–1336.
KAWASE, H., SATOH, T., FUKUTAKE, K. and 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.
LEE, M. K. W. and FINN, W. D. L., DESRA-2, Dynamic Eﬀective Stress Response Analysis of Soil Deposits
with Energy Transmitting Boundary Including Assessment of Liquefaction Potential (Soil Mechanics, 38,
Department of Civil Engineering, Univ. of British Columbia, Vancouver 1978).
LOMNITZ, C., CHAVEZ-GARCIA, F. J. and BARD, P. -Y. (1995), Two Articles on Gravity Waves; Discussion
and Reply, Bull. Seismol. Soc. Am. 85, 4, 1271–1274.
MARMARELIS, P. Z. and MARMARELIS, V. Z., Analysis of Physiological Systems, The White-Noise Approach(Plenum Press, New York and London 1978).
O’CONNEL, D. R. H. (1999), Replication of Apparent Nonlinear Seismic Response with Linear Wave
Propagation Models, Science 283, 5410, 2045–2050.
SATO, K., KOKUSHO, T., MATSUMOTO, M. and YAMADA, E. (1996), Nonlinear Seismic Response and Soil
Property during 1995 Hyogoken Nanbu Earthquake, Soils and Foundations, Special Issue on
Geotechnical Aspects of the January 17, 1995 Hyogoken Nanbu earthquake, 41–52.
SEED, H. B., WONG, R. T., IDRISS, I. M. and TOKIMATSU, K., Moduli and Damping Factors for Dynamic
Analyses of Cohesionless Soils (Report no. UCB/EERC-84/14, September 1984).
SOEDA, Y., TAMAI, H., NAKATSU, N., TAKEZAWA, K. and MAEGAWA, F., Nonlinearity and Irregularity on
Strong Seismic Motions by Vertical Arrays during the 1995 Kobe Earthquake. In Earthquake GeotechnicalEngineering (ed. Seco e Pinto) (Balkema, Rotterdam 1999) pp. 225–230.
SUN, J. I., GOLESORKHI, R. and SEED, H. B., Dynamic Moduli and Damping Ratios for Cohesive Soils(Report no. UCB/EERC-88/15, August 1988).
ZVOLINSKII, N. V. (1982), Wave Processes in Non-elastic Media, Problems of Engineering Seismology 23,
4–19.
(Received September 18, 2001, accepted May 2, 2002)

To access this journal online:
http://www.birkhauser.ch