Bulletin of the Seismological Society of America, Vol. 92, No. 5, pp. 1952–1969, June 2002

Changes in Shear Moduli of Liquefied and Nonliquefied Soils during the
1995 Kobe Earthquake and Its Aftershocks at Three Vertical-Array Sites
by Olga Pavlenko and Kojiro Irikura

Abstract Acceleration records of the 1995 Kobe earthquake and its first aftershocks, provided by vertical arrays at Port Island (PI), SGK, and TKS sites, wereused to estimate changes in the shear moduli of subsurface soils caused by strong
ground motion. Shear-moduli estimates are based on stress–strain relations determined for consecutive time intervals during the mainshock and aftershocks. We
found that in the upper 0–13 m at PI (reclaimed fill above the groundwater level),the shear modulus decreased by ⬃80%–90% of its initial value owing to liquefaction.
Signs of its recovery appeared 3–5 min after the decrease in the intensity of strong
motion. At depths of 13–18 m (reclaimed fill below the water level) and 27–32.5 m(water-saturated sand) at PI, 0–11 m at SGK, and 0–14 m at TKS (water-saturated
alluvium), the shear modulus decrease was estimated at about ⬃80%, 50%–60%,
60%–70%, and ⬃50%, respectively. Recovery began immediately after the decreasein the intensity of strong motion and was completed within a few minutes. The results
show that the behavior of the upper layers at sites PI and SGK, where the highest
accelerations were recorded, was dynamically unstable. So far as the method can beapplied to any site where vertical-array records and profiling data are available, accumulation of representative data on nonlinear soil behavior during strong motion
for various soil conditions is desirable. Summarizing and generalizing these data, wemay be able to predict soil behavior in future earthquakes.
Introduction
Although studies concerning the dynamic properties ofsoils date back several decades, they came to a qualitatively
new experimental level during the 1960s–1970s after thecatastrophic earthquakes in 1964 in Niigata (Japan) andAnchorage (Alaska). These earthquakes were most destructive and accompanied by extensive liquefaction of watersaturated, sandy soils. Another strong earthquake that occurred in 1971 in California greatly damaged the San
Fernando dam. After these events, dynamic testing of soils
became of priority importance, and a breakthrough wasmade by Professor Seed of Berkeley University and his
group (Seed and Lee, 1966). Their work initiated a wide
study of liquefaction phenomena in sandy soils (Casagrande,
1975; Castro, 1975; Ishihara et al., 1975). For reliable calculation of the vibration of structures based on soil layers,
dynamic (elastic, filtering, and damping) characteristics ofsoils were studied, and dynamic stability of soils was estimated. In particular, changes in their properties under dynamic loading were investigated, such as elastic-moduli degradation, increase in deformability, and decrease in rigidity.
The problem of the estimation of dynamic instability ofsoils (i.e., a higher probability of the deformability increase
and shear-modulus decrease in soils under dynamic, rather

than static, loading [Voznesensky, 1994]) is urgent. The examples of recent, destructive earthquakes (Spitak, 1988;
Loma Prieta, 1989; Northridge, 1994; Kobe, 1995; Neftegorsk, 1995; Taiwan, 1999; and Gujarat, 2001) show thatdamage to buildings and structures is high, and thus, morereliable and adequate estimates of stability of ground basements under dynamic loading are required.
Among recent methods of soil testing, laboratory dynamic methods are most diverse and numerous. The main
problem in laboratory testing is to create adequate dynamic
loading on a soil sample with the help of a laboratory installation. However, it seems to be impossible in principle
to create adequate loading (Lee, 1976). Among field methods, the geotechnical field methods, which aim at an estimation of soil rigidity under conditions of dynamic loading,such as penetration methods, dynamic sounding, pile testing,and so forth, seem to be more promising compared with theseismoacoustic and vibration methods. Great success has
been achieved in this field, and extensive data have been
accumulated concerning the behavior of cohesive and cohesionless soils under loadings of various intensities.
Nonetheless, disagreement remains between the in situ
conditions of dynamic loading on soils in earthquakes and1952

Changes in Shear Moduli of Liquefied and Nonliquefied Soils during the 1995 Kobe Earthquake

their laboratory and field-experimental substitutes. Soils
have been found to be such complex, natural, multicomponent systems that not only do they change the parameters ofpropagating seismic waves but also they are changed themselves. These phenomena are inseparable; they are rather
complicated and connected with the movement of underground water. Therefore, any data on soil behavior duringstrong ground motion and on the elastic moduli of soils obtained in situ are of particular importance. In situ variationsof the shear modulus from one observational point to anothercan be evaluated by the standard spectral ratio and horizontalto-vertical spectral ratio techniques (Beresnev and Wen,
1996; Dimitriu et al., 2000, 2001). However, since the spectra are calculated over long intervals of accelerograms, theresults are smoothed and do not allow determination of theshear-modulus variations during strong motion. Such information can be obtained from vertical-array records of strong
motion (Pavlenko and Irikura, 2001a, b). Seismic verticalarrays provide valuable data on soil behavior at different
depths during strong ground motion.
Vertical-array records of the 1995 Kobe earthquake undoubtedly testified to the soil nonlinearity at many recordingsites (Sato et al., 1996; Yasuda et al., 1996; Aguirre andIrikura, 1997; Ishihara et al., 1998; Kazama et al., 1998;
Kokusho and Matsumoto, 1998; Soeda et al., 1999; Pavlenko and Irikura, 2001a, b). Methods of backcalculation,
multireflection analysis, and inversion of records of the 1995
Kobe earthquake at Port Island (PI), SGK, TKS, and KNK
vertical-array sites were used to estimate the optimum valuesof S-wave velocity (Vs) and damping ratio in sedimentarylayers during the mainshock and an aftershock (Sato et al.,1996). The optimized values of Vs during the mainshock
were found to be smaller than their initial values obtainedby the PS-wave logging, especially for the upper liquefied
layer at PI. To trace the development of liquefaction, therecords at PI were divided into two parts, and it was foundthat liquefaction was still developing after the first period,
4.1–9.22 sec. Analysis of an aftershock with a maximum
surface acceleration of ⬃10 gal (cm/sec2) at PI, which occurred 2 min after the mainshock, showed that Vs for Pleistocene and Holocene layers underlying the reclaimed fill had
returned to their initial values, while Vs in the fill layer had
almost the same reduced values as in the mainshock, indicating that the fill was still under a liquefied condition (Satoet al., 1996). However, the computed horizontal acceleration
was slightly different from the observed values and quite
different at the level of GL-25 m for site TKS. Hence, it wasdifficult to simulate nonlinear seismic amplification by using
multireflection analysis with optimized properties for a longduration time history, including completely different types
of waves at site TKS (Sato et al., 1996).
The work was continued by Kokusho and Matsumoto(1998). Maximum acceleration and maximum shear strain
were estimated as a function of depth. A clear strain dependency was found in the backcalculated shear modulus and1953

damping, which is consistent with the laboratory test resultsfor each kind of considered soil. It was concluded that further study was needed for a better estimation of nonlinear
response during strong earthquakes.
The seismic behavior of an alluvial clay layer underlying reclaimed soil at PI was studied by Kazama et al. (1998).
The shear stress–strain relationship of the clay layer was
estimated from the mainshock by calculating separately theshear-strain and shear-stress evolution from filtered (f ⱕ 1.2
Hz) acceleration records above (GL-16 m) and below (GL32 m) the clay layer. Stress–strain relations corresponding
to low-frequency components were estimated for six successive time intervals within ⬃10 sec of strong motion.
However, strong-motion records by vertical arrays provide us with all the necessary information to reconstruct
stress–strain relations in situ, from the surface down to thelocation of the deepest recording device (Pavlenko and Irikura, 2002a). The obtained stress–strain relations describe
in detail soil behavior during strong motion and allow further
study of the qualitative and quantitative characteristics ofelastic nonlinearity of soils in different layers (Pavlenko and
Irikura, 2002b).
To determine the stress–strain relations in layers fromthe surface down to the location of the deepest device, thepropagation of a vertically incident seismic S wave in thehorizontally layered sediments is calculated. Assumptions
based on the profiling data are made about types of stress–
strain relations in different layers, and sets of various stress–
strain curves are generated. Item-by item examination is thenapplied to find groups of curves that show the best-fit approximation between the observed and simulated records.
As an example, records of the 1995 Kobe earthquake
provided by vertical arrays at PI, SGK, and TKS sites wereused to estimate the stress–strain relations in soils. A good
agreement was obtained between the observed and simulated
records, and a better understanding of processes occurringin soils during strong motion was achieved (Pavlenko and
Irikura, 2001a). It was shown that during the mainshock, themost significant (among the three studied sites) temporalchanges in soil behavior took place at sites PI and SGK,
where the strongest accelerations were recorded (⬃600 gal
at ⳮ83 m and ⬃340 gal on the surface at PI, and ⬃300 gal
at 97 m and ⬃650 gal on the surface at SGK). The obtained
stress–strain relations were found to vary with time at thesesites in surface layers (0 to 10–14 m) but remained stable indeeper parts. At PI, a progressive reduction of the shear modulus in the upper layers (0–13 m) due to liquefaction can betraced by stress–strain relations changing with time. At the
TKS site, the acceleration level was lower (less than ⬃100
gal at 100 m and ⬃200 gal on the surface). Soil behavior in
the upper layers can be described by a single, nonlinear
stress–strain relation, and changes in the soil behavior are
due to changes in the working areas of this curve, which
depend on the intensity of the imposed seismic motion. Since
soil conditions are similar at the three sites, the conclusion

1954

O. Pavlenko and K. Irikura

was made that the more complex soil behavior at PI and SGK
was due to the more intense input motion at these sites (Pavlenko and Irikura, 2002a).
The 1995 Kobe earthquake was accompanied by several
aftershocks just after the mainshock. In the present work, thedescribed method is applied to estimate stress–strain relations during weaker motion, that is, in the aftershocks occurring within 5 min after the mainshock. We used the obtained stress–strain relations for aftershocks to estimate
shear moduli at different depths at PI, SGK, and TKS sites.
Temporal changes in shear moduli during the mainshock andsubsequent aftershocks were traced, and estimates of a decrease in shear modulus are presented for PI, SGK, and TKS
sites.
It should be mentioned that aftershocks were superimposed on the background of basin-induced surface waves,which occurred after the mainshock. Basin-induced waves
recorded at PI and SGK sites after the mainshock can betreated as dynamic loading, and therefore, shear-modulus
estimates for aftershocks were obtained in a dynamic regime.

Data and Methods
Records of the mainshock and the subsequent aftershocks recorded at sites PI, SGK, and TKS were studied. Allof the sites are located around Osaka Bay and are close to
the epicenter of the 1995 Kobe earthquake. The distancesfrom PI, SGK, and TKS recording sites to the fault plane are
2 km, 6 km, and 16 km, respectively. Figure 1 shows thesoil stratigraphy and the arrangement of recording devices
at these sites. P-wave velocity indicates the groundwaterlevel at ⬃13 m at PI and 3–5 m at SGK and TKS sites.
The method of estimation of nonlinear stress–strain relations in layers from the surface down to the location of the
deepest device has been described by Pavlenko and Irikura(2002a). To determine the stress–strain relations in different
layers during strong ground motion, the medium from thesurface down to the location of the deepest device wasdivided into groups of layers, for which various types ofstress–strain relations were assumed according to the profiling data (Fig. 1). Three main types of stress–strain curves
were considered: (1) curves similar to those obtained in laboratory experiments by Hardin and Drnevich (1972) to describe the behavior of dense soils at depth (18–27 m and 33–
83 m at PI; below 11 m at SGK; and below 14 m at TKS);(2) curves of the “soft” type (Zvolinsky, 1982), similar totype 1 but with a much steeper slope, which are characteristic of liquefied soils (upper 0–13 m above the level ofunderground water at PI); (3) curves of the “hard” type, declining to the stress axis at large strains which are characteristic of water-saturated, sandy soils (Zvolinsky, 1982)(13–18 m and 27–33 m at PI; upper 0–11 m at SGK; upper
0–14 m at TKS). Thus, the choice of curves was physicallyjustified. To account for temporal changes in the soil behavior, records were divided into intervals of 1.5-sec duration.

Within each interval, the stress–strain relations were assumed to be stationary and to vary for different intervals.
Calculations were performed successively, interval by interval, and groups of best-fit stress–strain curves were foundfor each interval for two horizontal components at PI, SGK,
and TKS sites. The calculation algorithm based on the
lumped-mass method (Joyner and Chen, 1975) was used.
In this work, the method was applied to aftershocks ofthe 1995 Kobe earthquake that occurred within 5 min after
the mainshock, that is, to weak-motion records. Weak signals were found more difficult to analyze than strong ones,
because in weak motion, the working areas of the stress–
strain curves are located close to zero. As a result, many
curves look similar, and the difference between records simulated with different curves is very small. Therefore, we usedlarger data lengths than in strong motion (10–40 sec insteadof 1.5 sec), and we obtained the best-fit stress–strain curves
averaged over these time intervals. However, parameters of
the neighboring curves were close to each other, and inaccuracies in their estimation related to inadequate choice of
curves were smaller than in the case of the mainshock.
Profiling data used for nonlinear simulation at site PI
were taken from Aguirre and Irikura (1997). Data on shearwave velocity, shear-modulus degradation, and maximum
shear strain at all depths at SGK and TKS sites were takenfrom Soeda et al. (1999) and were applied to calculate theshear stress in failure, smax, and the low-strain shear modulus, Gmax, in different layers, according to Seed et al. (1984)
and Sun et al. (1988).
The obtained stress–strain relations in the mainshock
and aftershocks were used to estimate shear moduli at different depths. We calculated Gsec, which approximates average values of Gtan over the entire stress–strain loops (Kramer, 1996). As explained previously, one stress–strain curve
was used to describe the behavior of soils within a certain
group of layers at a given time interval. However, the working areas of stress–strain curves differ from layer to layer
within a group (as shown in Figs. 3b, 4b, and 5b), and estimates of Gsec vary within each group from layer to layer.
Weighted-mean values of Gsec were determined, which characterize average values of the shear modulus within a given
group of layers at a given time interval.
Figure 2 shows accelerograms of the mainshock together with the first aftershocks at PI, SGK, and TKS sites.
Long-period waves recorded in horizontal components at PI
and SGK sites may be standing waves due to S-wave reverberation in the basin sediments, probably including basininduced surface waves (Hatayama et al., 1995; Pitarka et al.,1998). Both sites PI and SGK belong to the Osaka basin,
and the corresponding depths to the basement are 1,000–
1,200 m and 1,200–1,400 m, respectively (Hatayama et al.,1995). Amplitudes of these waves achieved ⬃5 cm/sec2 atsite SGK and ⬃3 cm/sec2 at PI, and the period was about 4–
6 sec. Basin-induced waves were not recorded at site TKS
because of its location outside the limits of the Osaka basin,on the other side of the fault plane (Hatayama et al., 1995).

1955

2

VP

0

200

4.0

VP

8.0

2000

600

2

VS , m/s
2
VP, m/s

SGK-site

Failure strength, dyn/cm 2

VS

1000

400

Seismograph (GL 0,-24.9,-97.0 m)

Clay100. 0
0.0.0
0

80. 0

Clay

Clay60. 0

Clay40. 0

20. 0

0. 0

0

40 0

VP

8.0

2

2

m/s
S
2000 VP, m/

VS, m/m/sS
60 0

TKS-site

Failure strength, dyn/cm 2

4.0

1000

VS

20 0

Seismograph (GL 0,-25,-100 m)

0.0

0

0

Figure 1. The ground profiles and characteristic types of stress–strain relations at different depths at sites PI, SGK, and TKS. The
profiling data for PI were taken from the article by Aguirre and Irikura (1997). Data on shear-wave velocity, shear-modulus degradation,and maximum shear strain at sites SGK and TKS were taken from the article by Soeda et al. (1999); these data were used to calculate
shear stress in failure, smax, and the low-strain shear modulus, Gmax, in layers.

Port Island site2

VS , m/s
21500 V , m/s
P
600

1E+1 1E+2 1E+3 1E+4

VS

1000

400

Failure strength, dyn/cm
Seismograph (GL 0,-16,-32,-83 m)

1

500

0

0.1

200

0

1956

5:47:47

5:48:47

5:51:47EW-3.0
5:46:47

-1.5

0.0

1. 5

3. 0

a, m/s2-3.0
5:46:47

-1.5

0.0

1. 5

5:47:475:47:47

5:49:4 7

5:49:4 7

5 :50:47

5 :50:47

5:51:47

5:51:47

Time

5:52:47

UD

5:52:47-1.5
5:46:520.0

1.5a, m/s2-1.5
5:46:520.0

5:47:52

5:47:52

5:47:52

5:48:52

5:48:52

5:48:52

5:49:52

5:49:52

5:49:52

SGK

5:50:52

5:50:52

5:50:52

Time

5:51:52

UD

5:51:52

NS

5:51:52EW-1.0
5:46:510.0

1.0a, m/s2-1.0
5:46:510.0

1.0a, m/s2-1.0
5:46:510.0

1.0a, m/s2

5:47:51

5:47:51

5:47:51

5:48:51

5:48:51

5:48:51

5:49:51

5:49:51

5:49:51

TKS

Figure 2. Acceleration records of the 1995 Kobe earthquake and its aftershocks at sites Port Island (PI), SGK, and TKS at depthsof 32 m, 24.7 m, and 25 m, respectively. Reverberated S waves and basin-induced surface waves are seen in horizontal components
at sites PI and SGK.

5:48:47

5:48:47NS-1.5
5:46:521.5

5:52:47

3. 0

5 :50:47a, m/s2

5:49:4 7a, m/s2-3.0
5:46:47

-1.5

0.0

0.0

1.5

3. 0

1. 5

a, m/s2

a, m/s2

Port Island

5:50:51

5:50:51

5:50:51Time

5:51:51UD

5:51:51NS

5:51:51

EW

Changes in Shear Moduli of Liquefied and Nonliquefied Soils during the 1995 Kobe Earthquake

Results and Discussion
At PI, the strongest horizontal accelerations were measured during the mainshock (⬃500 gal for E–W and ⬃600
gal for N–S at 83 m), and liquefaction occurred in the upper
layers (0–13 m). First aftershocks were recorded immediately after the mainshock, when the surface soils were stillin a liquefied state. Figure 3a shows the observed and simulated records of the first five aftershocks at PI. Figure 3b
represents the corresponding stress–strain relations to depthsof 43 m, which characterize soil behavior during the mainshock and aftershocks. In the case of the mainshock, thestress–strain relations were obtained for 10 successive 1.5sec time intervals, whereas for the aftershocks, each groupof curves corresponds to one aftershock. The whole set ofobserved and simulated records of the mainshock at the threesites is presented in Pavlenko and Irikura (2002a).
Progressive liquefaction within the upper 0–13 m during
the mainshock can be traced by the evolution of stress–strain
relations: their slope decreases with increasing time, indicating a substantial reduction in shear modulus to approximately near-zero values (Fig. 3b). Stress–strain relations obtained for aftershocks at PI also show large slopes for theupper 0–13 m; however, these slopes are slightly steeper
than those obtained for final parts of the mainshock. Thiscan indicate the beginning of shear-modulus recovery a few
minutes after the mainshock.
At sites SGK and TKS, the imposed accelerations were
not as high as at PI (⬃300 gal and ⬃100 gal, respectively).
Figure 4a and 4b shows the observed and simulated aftershocks and the corresponding stress–strain relations, respectively, for the mainshock and aftershocks at site SGK.
Changes in soil behavior occurred in the upper 0–11 m andcan be traced by changes in the shapes of the stress–strain
curves. Stress–strain relations characterizing the soil behavior in the aftershocks show a complete recovery of the shear
modulus for the upper layers. At site TKS, the soil behaviorof the upper 0–14 m during the mainshock and aftershocks
can be described by a single stress–strain relation. Changes
in the shear modulus are due to changes in the working areasof this curve, depending on the intensity of the imposed motion. The observed and simulated accelerograms of aftershocks at site TKS are presented in Figure 5a, and the corresponding stress–strain curves for the mainshock and theaftershocks are shown in Figure 5b.
The difference between aftershocks recorded in liquefied (PI) and nonliquefied (SGK and TKS sites) soils isclearly seen in Figures 3a, 4a, and 5a. At sites SGK and
TKS, acceleration amplitudes on the surface were noticeably
higher than at the PI site. A “natural” amplification of seismic waves in the upper layers was observed at SGK and
TKS sites (Figs. 4a and 5a), whereas at PI, where soils wereliquefied, the acceleration amplitudes were very small on thesurface for all aftershocks except the last one, which showed
higher amplitudes. This confirms that 5 min after the ter-

1957

mination of strong motion, recovery of shear moduli in liquefied layers was already noticeable.
Figure 6a–6e shows summary plots of shear moduli inlayers 0–13 m (a), 13–18 m (b), and 27–32.5 m (c) at PI; 0–
11 m at SGK (d); and 0–14 m at TKS (e) as estimated fortwo horizontal components during the mainshock and aftershocks. These estimates are based on the corresponding
stress–strain relations shown in Figures 3b, 4b, and 5b. Below the shear-moduli curves, the corresponding observed
and simulated acceleration-time histories of the mainshock
are presented. Simulated records correspond to shear-moduli
evolution. A good agreement between the observed and simulated records testifies to the validity of the obtained shearmoduli estimates.
To evaluate the accuracy of the obtained shear-moduli
values, we analyzed two main sources of possible errors: (1)errors introduced by inexact determination of stress–strain
relations in the layers and (2) scattering of shear-modulus
estimates within groups of layers. Errors introduced by inexact determination of stress–strain relations in layers canbe roughly evaluated by deviations of the simulated records
from the observed ones. For the mainshock, these meansquare deviations, normalized by the intensities of the observed records, are ⬃0.3–0.4 for PI for frequencies up to
12 Hz and for site SGK, for frequencies up to 20 Hz, andslightly higher, ⬃0.5–0.6 for TKS site, for frequencies up
to 20 Hz. For aftershocks at PI, SGK, and TKS, these values
are higher: ⬃0.8–1.3, ⬃0.6–0.8, and ⬃0.7–1.2, respectively,
for frequencies up to 20 Hz. These values and accelerationtime histories show a rather good agreement between theobserved and simulated motions. Thus, although the errors
introduced by inexact determination of the stress–strain relations in layers can reach tens of percent, accounting forthe large numbers of points and the conformity of estimatesobtained for E–W and N–S components, we conclude thatthe main trends of shear-moduli changes were found correctly. Scattering (i.e., variations) of shear-modulus estimates within groups of layers is due to the fact that although
the stress–strain relation is the same for all of the layers
within the group, the working intervals are different for different layers, depending on elastic parameters of the layers.
Therefore, shear-modulus estimates vary within each group,
and we used the mean values and calculated confidence limits, which are shown in Figure 6a. They are not shown inother figures because of their smallness.
PI SiteIn the upper 0–13 m at PI, shear modulus decreased
gradually during the mainshock (Fig 6a). The decrease wasmost noticeable in the first 3–5 sec of strong motion, whenit achieved ⬃ 80%–90% of the initial value. Shear modulus
reached its minimum during the following 5–7 sec andshowed a tendency to increase a few minutes after a decreasein the intensity of the strong motion. Scattering of shearmodulus estimates corresponding to aftershocks indicatedthat recovery of shear modulus proceeded nonuniformly,

1958

O. Pavlenko and K. Irikura

Figure 3. (Caption on next page.)Changes in Shear Moduli of Liquefied and Nonliquefied Soils during the 1995 Kobe Earthquake

1959

Figure 3.

Acceleration-time histories of aftershocks in PI, observed and simulated(a), and the obtained stress–strain relations for the mainshock and aftershocks at PI (b).
The stress–strain curves are scaled in relative units, which are the same for all timeintervals at a given depth. Stress and strain are scaled in the manner described by Joyner
and Chen (1975); stress was normalized by multiplying by 1/smax, and strain was normalized by multiplying by Gmax /smax .
(Continued on next page.)which can be attributed to the movement of groundwater.
The upper subsoil layer is a decomposed granite soil called“Masado,” which is a sandy soil with much gravel, silt, andclay, and the contents of gravel and fines are ⬃0%–65% and5%–35%, respectively (Yasuda et al., 1996). According tomany researchers, the presence in soil of even a small quantity of clay particles (1.5%–2% by Gumensky, 1965) imparts
thixotropic properties to the soil (Boswell, 1949; Akkerman,
1958; Gumensky, 1965). Soils like Masado, which containsmall quantities of fines, possess intermediate properties between noncohesive and cohesive (quasi-thixotropic) soils.

For such weakly cohesive soils, a quick decrease in rigidity even for rather low levels of dynamic loading is typical, which is due to their mixed composition. In such soils,
the Coulomb friction is less than in a pure sand, and thecoagulation lattice is weak and usually discontinuous. The
decreased permeability hampers the dissipation of pore pressure and facilitates liquefaction. In such soils, gradual andslow hardening after termination of dynamic loading includes two processes, such as weak compression of the soil
accompanied by the loss of water and thixotropic strengthening of the coagulation lattice. The compression proceeds

1960
Figure 3. (Caption on previous page.)Changes in Shear Moduli of Liquefied and Nonliquefied Soils during the 1995 Kobe Earthquake

Figure 4. Acceleration-time histories of aftershocks at site SGK, observed and simulated (a), and the obtained stress–
strain relations for the mainshock and aftershocks at SGK (b). The axes scales of the stress–strain curves are in relative
units, which are the same for all time intervals at a given depth. Stress and strain are scaled as in Figure 3b.

1961

1962
Figure 4. (Caption on previous page.)Changes in Shear Moduli of Liquefied and Nonliquefied Soils during the 1995 Kobe Earthquake

Figure 5. (Caption on next page.)

1963

1964

O. Pavlenko and K. Irikura

Figure 5.

(Continued) Acceleration-time histories of aftershocks at site TKS, observed and simulated (a), and the obtained stress–strain relations for the mainshock
and aftershocks at TKS (b). The axes scales of the stress–strain curves are in relative units, which are the same for all time intervals at a given depth. Stress and strain
are scaled as in Figure 3b.

in accord with the dissipation of pore pressure, and its rateis determined by the filtering properties of a soil, whereas
the thixotropic strengthening takes several hours (maximum,
1 day) and proceeds rather nonuniformly during the first few
minutes owing to continued water movement (Voznesensky,1999). This recovery mechanism consisting of the two processes, pore pressure dissipation and thixotropic strengthening, seems to be probable at PI, considering the soil composition, the nonuniformity of recovery during the first few
minutes (Fig. 6a, 6b, 6d, 6e), and the total time for recovery,
which was estimated by Aguirre and Irikura (1997). They
showed by means of S-wave-velocity inversion and simulation of the S-wave propagation accounting for the soil non-linearity, that (1) the rigidity of the upper layers decreased
sharply during the largest ground shaking before liquefaction happened; (2) the biggest change in the soil rigidity wasobserved several seconds after the occurrence of liquefaction
and then it gradually recovered; (3) the most noticeablechanges in shear moduli occurred in the layer from 5 to 16
m depth; and (4) the aftereffects of liquefaction remained atleast 3 hours after the mainshock but no more than 24 hours.
Figure 6b shows the changes in shear modulus at depthsof 13–18 m and the corresponding observed and simulated
accelerograms at 16 m. Changes in the shear moduli in thislayer during the first few seconds of strong motion were aslarge as in the surface layer. The composition of soils at 13–

1965

Figure 5. (Caption on facing page.)

1966

O. Pavlenko and K. Irikura

Figure 6. Shear-moduli evolution of subsurface soils at sites PI, SGK, and TKS
during the 1995 Kobe earthquake and its aftershocks.

Changes in Shear Moduli of Liquefied and Nonliquefied Soils during the 1995 Kobe Earthquake

18 m was similar to that in the upper layers, but the soils
were saturated with water (Fig. 1). According to our computations, a decrease in rigidity in these layers started withthe beginning of strong motion and continued for 3–4 sec,that is, until the intensity of the imposed motion decreased.
After that, compression and probably thixotropic strengthening took place. The maximum decrease of shear modulus
was about 80%, and the obtained stress–strain relations show
its recovery during the following few minutes (Fig. 6b). The
same tendency of shear-modulus changes can be revealed
from stress–strain curves presented in Kazama et al. (1998),
that is, its decrease during the first 1–4 sec of strong motion
and the following recovery during the next 6–9 sec. If weaverage the obtained shear-modulus changes in 13–18 m,
our results also agree with the results obtained in Sato et al.
(1996), who identified an 80% and 50% reduction of Vs forthe saturated and unsaturated fill layers, respectively.
At depths of 27–32.5 m at PI, shear modulus decreased
by 50%–60% of the initial value during the first few secondsof strong motion (Fig. 6c), which agrees on the whole withestimates made by Sato et al. (1996) that the optimized Vs
values for the total duration of the mainshock (20 sec) were⬃20% and 40% smaller than the initial values obtained by
the PS-wave logging test for the Pleistocene (gravel/clay,
50–83 m) and Holocene (sand/clay, 20–50 m) layers, respectively.
SGK and TKS Sites
At site SGK, a shear modulus decrease was also observed in the upper layers, namely, ⬃ 60%–70% of the initial value; recovery was completed during the following
minutes. Figure 6d shows the temporal changes of shear
modulus in the upper 0–11 m at site SGK and the observed
and simulated accelerograms of the mainshock on the surface. A sharp decrease in the rigidity of the upper layers atsite SGK coincided with the beginning of strong motion, as
at PI. Recovery began just after a decrease in the intensityof seismic motion and continued for a few minutes, until theshear modulus attained its final value. Soil conditions at PI
and SGK were similar; therefore, the shear-modulus recovery at site SGK also included the soil compression accompanied by the loss of water and probably, the thixotropic
strengthening of the structural coagulation lattice.
The role of thixotropic strengthening in the shear-modulus recovery at PI, SGK, and TKS sites is not clear. From
one side, the recovery curves (Fig. 6a–6e) show some regularities of hardening, which are considered typical for
weakly cohesive soils, such as the beginning of the recovery
just after a decrease in the intensity of strong motion (Fig.
6b, 6c, 6d, 6e), the nonuniformity of recovery during thefirst few minutes (Fig. 6a, 6b, 6d, 6e), and the rate and total
time of the recovery. In the cases shown in Figure 6b, 6c,
6d, and 6e, recovery began just after a decrease in the intensity of strong motion. It is known that “instantaneous thixotropic strengthening” occurs just after termination of dynamic loading because the number of free particles taking

1967

part in the reconstruction of the structural lattice is maximum
at that time. The highest increase in rigidity is observed during the first 3–10 minutes, and the presence of water accelerates this process because it increases the mobility of particles and their ability to reorient and establish the mostenergetically advantageous contacts (Voznesensky, 1999).
Another peculiarity of some curves of the shear-modulus
recovery is a disagreement between the initial and final values of rigidity (Fig. 6c, 6d). Soils do not possess an important property of thixotropic systems, such as the entire reversibility of the loss and recovery of their rigidity. Their
final rigidity achieved after recovery is usually either less or
higher than the initial one, depending on the inhomogeneityof a soil and the intensity of the preceding dynamic loading.
Recovery of the rigidity in soils is related to phenomena
aimed at ordering the system by decreasing its entropy, suchas (1) increasing the energy of particle interaction, (2) destruction of structural bonds with low-level activation energy, (3) increasing the number of contacts between particles
within a unit of volume (Voznesensky, 1999), and (4) redistribution of water to make it more uniform (Barber,1958). These effects appear more clearly in cases of stronger
dynamic loading, only in the first cycle of “destruction–
recovery.”
However, from the other side, we have neither the results of laboratory experiments, proving that “Masado” doesbear thixotropic properties, nor any further evidence. The
nonuniformity of the shear-moduli recovery and the disagreement between the initial and final values of the shear
moduli can be the consequence of inaccuracies of theshear-moduli estimates. Therefore, in our case, thixotropic
strengthening is only probable, although we cannot discard
thixotropic mechanisms, because in weakly cohesive soils,
consolidation is usually accompanied by thixotropic
strengthening of the coagulation lattice, and the separationof compression and thixotropic effects is virtually impossible (Voznesensky, 1999).
Thus, our results confirm that the upper 10–13 m at PI
and SGK are dynamically unstable systems. At the sametime, the sensitivity of soils to dynamic loading obviously
decreases with depth, and no changes were observed indeeper layers at sites PI and SGK. According to Voznesensky (1999), the largest decrease in rigidity was observed at
the acceleration levels of 0.5–0.7g, and these values wereattained at both PI and SGK sites. On the other hand, at site
TKS, the highest measured accelerations were ⬃ 0.2g, andthe changes in soil rigidity of the upper layers were not sohigh. A single stress–strain relation was used to describe thesoil behavior in the upper 0–14 m at site TKS, and changes
in the shear modulus (⬃50% of the initial value) were due
to changes in the working areas of this stress–strain relation,
depending on the intensity of the imposed motion (Fig. 6e).
These changes occurred mostly in the upper meters below
the water level, at 4–7 m depth (Fig. 5b).
Our estimates of the shear-modulus decrease at sites
SGK and TKS differ from those obtained by Sato et al.

1968

O. Pavlenko and K. Irikura(1996). They found that at site SGK the Vs optimized for the
total duration of the mainshock was ⬃25% and 45% smaller
than the initial value for the Pleistocene and Holocene layers,respectively. At site TKS, sand boils were observed after thequake, indicating that liquefaction took place at this site. The
optimized Vs was ⬃10% and 30% smaller than the initialvalue for the Pleistocene and Holocene layers, respectively.
The optimized velocity of the sand layer at a depth of ⬃10m was reduced by 80% from the initial value, as in theliquefied-fill layer at PI (Sato et al., 1996). In contrast to
their results, ours show stability in the shear moduli at depths
below 11–14 m at sites SGK and TKS and no liquefaction
at site TKS, yet liquefaction at TKS would be very probablefor larger strains. Because the accuracy of the simulation in
Sato et al. (1996) is not very high, as they themselves admit,
and because their simulated records differ substantially fromthe observed ones, especially at a depth of 25 m, we consider
our estimates to be more exact.
Shear-moduli estimates at sites PI and SGK were obtained in a dynamic regime because of the presence of reverberated S waves and basin-induced surface waves. Thesewaves, representing standing waves in the Osaka basin, arenoticeable at all depths at sites PI and SGK (Fig. 3a and 4a)and possess amplitudes of 0.03–0.06 m/sec2 at a predominant frequency of ⬃0.2 Hz. The loading on soils created bythese waves should be treated as dynamic, rather than quasistatic, because soil conditions are mostly undrained (Ishihara, 1998) when the soil permeability is not high enoughfor a complete dissipation of excess pore pressure during thephase of loading in each cycle. Dynamic loading on soils
created by these waves is not large when compared with themainshock, and most probably under this loading, soils cameto a dynamic equilibrium characterized by a closed hysteretic loop.

Conclusions
The method of estimating nonlinear time-dependentstress–strain relations from vertical-array records, which has
been described in Pavlenko and Irikura (2002a), was applied
to study the changes of shear moduli in different layers atsites PI, SGK, and TKS during the 1995 Kobe earthquake
and its first aftershocks.
Analysis of the aftershocks that occurred shortly after
the mainshock, as well as of the mainshock itself, showedthat the most noticeable changes in soil rigidity occurred in
the upper layers at PI. The maximum decrease was observed
to be ⬃80%–90% in the upper reclaimed-fill layer (0–18 m)
and ⬃50%–60% in water-saturated sands at 27–32.5 m.
However, the shear-modulus recovery began 3–5 min after
the decrease in the intensity of strong motion in the upper
reclaimed layers above the level of underground water (0–13 m) and just after the decrease in the intensity of strong
motion in the reclaimed fill below water level (13–18 m) andat depths of 27–32.5 m.
At sites SGK and TKS, changes in the shear moduli

occurred in the upper water-saturated alluvial and reclaimedsoil layers, 0–11 m and 0–14 m, and reached ⬃60%–70%
and ⬃50%, respectively. At site SGK, shear-modulus recovery began after the decrease in intensity of the imposed
motion and was completed in a few minutes. At site TKS,
where recorded accelerations were substantially lower thanat PI and SGK, no changes in the stress–strain relations were
noticed, in contrast to sites PI and SGK. A single stress–
strain relation can be used to describe the soil behavior in
the upper layers at site TKS and the changes in soil rigidity
were due to changes in the working areas of this curve, depending on the intensity of the imposed seismic motion.
Thus, soils in the upper layers at sites PI and SGK andsandy soils at depth at PI showed a dynamically unstablebehavior during strong motion; that is, they changed theirproperties by decreasing their rigidity and increasing deformability. Soil behavior during the 1995 Kobe earthquake
is a good illustration of the fact that soils are such natural
multicomponent systems that they not only change the dynamic characteristics of propagating seismic waves but thatthey also change their own properties, and these processes
are inseparable from each other.
The method of estimation of the stress–strain relations
from vertical-array data show in detail soil behavior duringstrong motions and changes in the shear moduli in different
layers during the mainshock and aftershocks. Such in situ
observations represent a valuable complement to laboratory
and field studies of soil properties.
Insofar as the method can be applied at any site wherevertical-array records and profiling data are available, accumulation of representative data on nonlinear soil behavior
under strong motion in different soil conditions is desirable.
Summarizing and generalizing these data, we may become
able to find the rational basis of how the model (the stress–
strain curves) would change if the input changes in one or
another site, that is, to predict soil behavior in future earthquakes.

Acknowledgments
We thank Yoshio Soeda, Hideki Tamai, and Naoyuki Nakatsu (Kansai
Electric Power Company) for acceleration records of the vertical arrays at
SGK and TKS recording sites. The acceleration records of the vertical arrayin PI were kindly provided by Kobe City. This study was partially supportedby Special Coordination Funds of the Science and Technology Agency of
the Japanese Government (PI: Kojiro Irikura).

References
Aguirre, J., and K. Irikura (1997). Nonlinearity, liquefaction, and velocity
variation of soft soil Layers in PI, Kobe, during the Hyogo-ken Nanbu
Earthquake, Bull. Seism. Soc. Am., 87, 1244–1258.
Akkerman, E. (1958). Thixotropy and fluidity of fine-grained soils, Prob.
Eng. Geol., 1, 73–92.
Barber, E. S. (1958). Discussion on the paper ‘Thixotropic characteristics
of compacted clays’ by Seed, H. B., and Chan, C. K. J. Soil Mech.
Found. Div.

Changes in Shear Moduli of Liquefied and Nonliquefied Soils during the 1995 Kobe Earthquake
Beresnev, I. A., and K.-L. Wen (1996). Nonlinear soil response—a reality?,
Bull. Seism. Soc. Am., 86, 1964–1978.
Boswell, P. G. H. (1949). A preliminary examination of the thixotropy ofsome sedimentary rocks, Q. J. Geol. Soc. London, 104, 499–526.
Casagrande, A. (1975). Liquefaction and cyclic deformation of sands: acritical review (lecture), 5th Panamerican Conference on Soil Mechanics and Foundation Engineering, Buenos Aires, 80–133.
Castro, G. (1975). Liquefaction and cyclic mobility of saturated sands,
Proc. ASCE, 101, GT6, 551–569.
Dimitriu, P., N. Theodulidis, and P. Y. Bard (2000). Evidence of nonlinear
site response in HVSR from SMART1 (Taiwan) data, Soil Dynam.
Earthquake Eng. 20, 155–165.
Dimitriu, P., N. Theodulidis, P. Hatzidimitriou, and A. Anastasiadis (2001).
Sediment non-linearity and attenuation of seismic waves: a study ofaccelerogramms from Lefkas, western Greece, Soil Dynam. Earthquake Eng. 21, 63–74.
Gumensky, B. M. (1965). Basis of physics and chemistry of clay soils andtheir use in civil engineering, Stroiizdat, Leningrad, Moscow, 255.
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.
Hatayama, K., K., Matsunami, T. Iwata, and K. Irikura (1995). Basininduced Love waves in the eastern part of the Osaka Basin, J. Phys.
Earth, 43, 131–155.
Ishihara, K., F. Tatsuoka and S. Yasuda (1975). Undrained deformationand liquefaction of sand under cyclic stress, Soils Found. 15, 29–44.
Ishihara, K., M. Cubrinovski, and N. Takeshi (1998). Characterization ofundrained behaviour of soils in the reclaimed area of Kobe, Soils
Found. (special issue on geotechnical aspects of the January 17 1995
Hyogoken-Nambu earthquake), 2, 33–46.
Joyner, W. B., and T. F. Chen (1975). Calculation of nonlinear groundresponse in earthquakes, Bull. Seism. Soc. Am., 65, 1315–1336.
Kazama, M., A. Yamaguchi, and E. Yanagisawa (1998). Seismic behaviorof an underlying alluvial clay on man-made islands during the Hyogoken-Nanbu earthquake, Soils Found. (special issue on geotechnicalaspects of the January 17 1995 Hyogoken-Nambu earthquake), 2,
23–32.
Kokusho, T., and M. Matsumoto (1998). Nonlinearity in site amplificationand soil properties during the 1995 Hyogoken-Nambu earthquake,
Soils Found. (special issue on geotechnical aspects of the January 17
1995 Hyogoken-Nambu earthquake), 2, 1–9.
Kramer, S. L. (1996). Geotechnical Earthquake Engineering, Upper Saddle
River, Prentice Hall, NJ, 230–238.
Lee, K. L. (1976). Influence of end restraint in cyclic triaxial tests. WES
CR S-76-1, USAE Waterways Experiment Station, Corps of Engineers, Vicksburg, MS.
Pavlenko, O. V., and K. Irikura (2002a). Estimation of nonlinear time-1969

dependent soil behavior in strong ground motions based on verticalarray data, Pageoph (in press).
Pavlenko, O. V., and K. Irikura (2001b). Types of elastic nonlinearity of
sedimentary soils, Geophys. Res. Lett., (in press).
Pitarka, A., K. Irikura, T. Iwata, and H. Sekiguchi (1998). Three-dimensional simulation of the near-fault ground motion for the 1995 Hyogoken Nanbu (Kobe), Japan, earthquake, Bull. Seism. Soc. Am., 88,
428–440.
Sato, K., T. Kokusho, M. Matsumoto, and E. Yamada (1996). Nonlinear
seismic response and soil property during 1995 Hyogoken Nanbu
earthquake, Soils Found. (special issue on geotechnical aspects of theJanuary 17, 1995 Hyogoken Nanbu earthquake), 2, 41–52.
Seed, H. B., and K. L. Lee (1966). Liquefaction of saturated sands during
cyclic loading, J. Soil Mech. Found. Div. 92, 105–134.
Seed, H. B., R. T. Wong, I. M. Idriss, and K. Tokimatsu (1984). Moduli
and damping factors for dynamic analyses of cohesionless soils, Report No. UCB/EERC-84/14, September.
Soeda, Y., H. Tamai, N. Nakatsu, K. Takezawa, and F. Maegawa (1999).
Nonlinearity and irregularity of strong seismic motions by verticalarrays during the 1995 Kobe earthquake, in Earthquake GeotechnicalEngineering, Seco e Pinto ed., Balkema, Rotterdam, 225–230.
Sun, J. I., R. Golesorkhi, and H. B. Seed (1988). Dynamic moduli anddamping ratios for cohesive soils, Report No. UCB/EERC-88/15,
August.
Voznesensky, E. A. (1994). Dynamic instability of soils and rocks, in Proceedings of the 7th International IAEG Congress, Lisboa, A. A. Balkema, Rotterdam, Brookfield, 683–692.
Voznesensky, E. A. (1999). Dynamic instability of soils, Editorial UPSS,
Moscow, 264.
Yasuda, S., K. Ishihara, K. Harada, and N. Shinkawa (1996). Effect of soil
improvement on ground subsidence due to liquefaction, Soils Found.
(special issue on geotechnical aspects of the January 17 1995 Hyogoken-Nambu earthquake), 2, 99–107.
Zvolinsky, N. V. (1982). Wave processes in non-elastic media, Probl. Eng.
Seism. 23, 4–19.
Institute of Physics of the Earth
Russian Academy of Sciences,
Moscow Russia(U.P.)
Disaster Prevention Research InstituteKyoto University,
Kyoto, Japan(K.T.)
Manuscript received 26 March 2001.