Bulletin of the Seismological Society of America, Vol. 87, No. 5, pp. 1244-1258, October 1997

Nonlinearity, Liquefaction, and Velocity Variation of Soft Soil Layers
in Port Island, Kobe, during the Hyogo-ken Nanbu Earthquakeb y Jorge Aguirre and Kojiro Irikura

Abstract Clear nonlinear behavior is analyzed from the acceleration records ofthe 1995 Hyogo-ken Nanbu earthquake at Port Island, Kobe. From four triaxial instruments placed at four different depths, the surficial effects during strong groundmotions were compared with those during weak motions before and after the mainshock. We used a spectral ratio technique and a nonlinear inversion for velocity
structure to analyze the data. From the spectral analysis, we observed a large variationof the spectral ratios between the surface and different depths during the strong
ground motions and during the liquefied state. The spectral ratios after the mainshock(i.e., after the liquefied state) are different from those before the mainshock. The
peak frequencies in the spectral ratios after the mainshock are shifted to lower frequencies with respect to those in the spectral ratios before the mainshock. We inverted the S-wave velocities using a genetic algorithm technique to determine thevelocity structure before, during, and after the mainshock. The S-wave velocity structure before and after the mainshock was found to be different. Specifically, the Swave velocity of the second layer (5 m to 16 m depth) after the mainshock was 20%
lower than before. Our analysis shows that the liquefied state remains at least 3 hr
after the mainshock but no more than 24 hr. The rigidity of the soil decreased close
to zero when liquefaction happened and later increases gradually following a trendthat resembles a consolidation curve. The strong influence of nonlinearity during themainshock yielded a big reduction of the horizontal surface ground motions, so thatthe observed horizontal peak acceleration was only about 25% of the peak acceleration expected from the linear theory. However, the nonlinear effects in the verticalpeak acceleration were not significant.

Introduction
Many authors have been studying the amplification produced from the sedimentary soils in order to predict the motions during large earthquakes, as well as to reproduce theground motions during large earthquakes at sites where theground motions have not been recorded. The amplificationof the seismic waves basically originates from the strong
contrast between the physical properties of the rocks and thesediments. To evaluate this amplification, conventionally,
the seismic response of the soils has been treated as linear
in the seismological community. That means strain andstress are related linearly by a constant rigidity independently of the level of strain. This assumption is acceptable
for low levels of strain (for 7c = 0.0001% and 0.001%,
where 7c is the cyclic shear strain amplitude, Vucetic andDobry, 1991), from soil laboratory testing results. But forlarger stress-strain levels, the soil laboratory testing results
showed a nonlinear relation, that is, the nonlinear characterof the soil response.
The basic reason that the seismologist conventionallyhas been using a linear relation is that for the level of strains
in most of the seismological observations the shear modulus
reduction is small. But, according to the laboratory test resuits, the nonlinear relation in the case of large strains may
play an important role in ground motions at soil sites near
the source during large earthquakes. Since the conditions in
the laboratory are not the same to those in the field, it isimportant to investigate the nonlinearity of the soil response
for large strain levels through in situ observations by separating the soil response from other factors such as sourceand path effects in seismic data.
Some authors have been trying to find observational
proof of the nonlinearity from seismological data and estimate how much it influences strong ground motions (Sugito
and Kameda, 1990; Chin and Aki, 1991; Ordaz and Faccioli,
1994; Aguirre et aI., 1994; Beresnev et al., 1995; Satoh etal., 1995). One of the problems that they have been facingis the lack of records for large strain levels during the strong
ground motions. Also, the mixed influence of the nonlinear

1244

Nonlinearity, Liquefaction, and Velocity Variation of Soft Soil Layers in Port lsland, Kobe
soil response to the recorded motion and other effects makes
it difficult to study accurately. Since the nonlinear effecttypically occurs in unconsolidated soft soils, our analysis
should be focused on very surficial layers. Then, in order tostudy the nonlinear effect from the seismological point ofview, it is necessary to extract the local site effects from theobserved records. There are several ways to overcome these
problems. The transfer function of seismic waves in soil layers is estimated by calculating the spectral ratio between the
motion at the soil surface and the motion at the basement
underlying the soil layers. In many studies, the rock outcrop
motion is used to estimate the motion at the basement. Thatis an approximation for the real motion at the basement under the soil layers; the accuracy of this approximation depends on the simplicity of the topographical conditions
around rock outcrop station. In most of the cases, the surface
topography at the rock outcrop site is very complex.
Another option for estimating the soil layers effect is touse data from a vertical array of seismometers. Using therecords from the vertical arrays, it is possible to separate thesite effect from the source and path effects and then verify
the behavior of the soil layers between the stations at different depths whenever an earthquake is recorded. So it ispossible to identify the nonlinear behavior and the changesof the physical properties in the layers with the lapsed time.
Unfortunately, the borehole stations are not so widely deployed as surface stations. This is the reason why the nonlinear behavior has been difficult to identify. Even if the
deeper borehole stations of the vertical array do not reachthe basement, it is useful to separate the local effects at sites
just where the nonlinear effects are expected to occur, in theunconsolidated soil layers closer to the surface.
Here, the Port Island vertical array data are used to study
the nonlinear response and the time variation of the velocities of the soil layers that may be caused by the liquefaction
during the intense part of the 1995 Hyogo-ken Nanbu earthquake as well as the variation of soil properties during theliquefaction state. First, the velocity structure model is precisely determined analyzing the five events before the mainshock by nonlinear inversion being revised from P-S waves
logging data. Next, the revised model is used to simulate themainshock for the linear and nonlinear cases. The long-duration record during the mainshock is used to discuss thevariation of the response within 6 rain after the mainshock.
Finally, the inversion of the velocity structure for the aftershocks was made. The purpose of this article is to discuss
the time variation of local site effects and its influence on
the ground motion.
A r r a y and D a t a
The array is composed of four triaxial accelerometers
placed at GL-0 m, GL-16 m, GL-32 m, and GL-83 m, where
GL means "ground level," hereafter identified as PR4, PR3,
PR2, and PR1, respectively. The original velocities, type ofsoil, and positions of the sensors are shown in Table 1 (De-1245

velopment Bureau, Kobe City, 1991). Since the horizontal
components of the instrument placed at 83 m depth deviate
19° in clockwise direction as studied by Ansary etal. (1995),
we correct all of the records that we used in this study, applying a rotation of 19 ° in the opposite direction.
The set of records used here is shown in the Table 2. It
consists of records from five earthquakes occurring in 1994
before the mainshock, the long-duration (6 rain) record of
the mainshock, and 12 aftershocks occurring 3 hr to 10
months after the mainshock. The epicentral locations of allof these earthquakes are shown in Figure 1. The record of
the mainshock with a duration of 360 sec was divided intonine sections, we refer to them as mainshock and sections
A to H. One example of this very long-duration record withthe sections A to H is shown in Figure 2. The records of the
N - S component at - 16 m show an unusual behavior during
the aftershocks from 1995/01/17 08:58 hr to 1995/05/19

Table 1
Original Soil Model in Port Island Vertical Array
Depth
(rn)

Soil T y p e

0-2.0
2.0-5.0
5.0-12.6
12.6-19.0
19.0-27.0
27.0-33.0
33.0-50.0

Location of
Accelerometers

Vp[km/sec]

V~[km/sec]*PR4-0 m

0.260
0.330
0.780
1.480
1.180
1.330
1.530

0.170

1.610

0.350
0.303
0.320

GravelR

SandygravelR
Clay
Sand
Sandygraveland sand50.0-451.0
Sand
6 1.0-79.0
Clay79.0-(85.0) Sandygravel*PR3-16m*PR2-32 rn*PR1-83m

2.000

0.210
0.180
0.245
0.305

R = Reclaimedland.

Table 2
Seismic Events Used in This Study
Date(year.month.day)Time(hr:rnin:sec)

Depth
(kin)

Magnitude

94.06.28
94.07.28
94.10.24
94.11.09
94.11.10
95.01.17
95.01.17
95.01.18
95.01.26
95.02.02
95.02.02
95.02.02
95.02.03
95.02.18
95.05.04
95.05.19
95.09.12
95.10.14

13:08:53.02
10:01:52.04
11:51:10.72
20:26:56.41
00:38:17.72
05:46:46.74
08:58:16.14
05:25:40.39
01:01:18.38
12:57:22.13
16:04:19.55
16:19:27.71
20:36:55.32
21:37:33.66
l 7:42:02.18
20:35:38.97
06:30:26.41
02:03:59.15

16.0
11.5
15.1
10.4
11.1
16.0
18.8
9.56
11.01
5.94
12.80
17.24
5.22
20.42
16.2
20.8
15.7
16.8

4.6
4.1
4.3
4.1
3.9
6.9
4.7
3.0
3.3
2.4
3.4
3.4
3.0
3.8
4.3
4.1
3.9
4.8

1246

J. Aguirre and K. Irikura134"30'135"00'135"30"136"00'35* 00"35* 00'34 ° 30'

34 ° 30'134"30'135 ° 00'

135 ° 30'

20:35 hr, probably because of a malfunction of the sensor.
This sensor was replaced by a new one around 20 May 1995.
Those irregular data will not be used in this analysis.
MethodWe analyze the variation of the S wave's velocity structure on Port Island vertical array before, during, and after
the mainshock by means of spectral ratios, correlation analysis, and nonlinear inversion of the velocity structure.
In the spectral ratio analysis, we choose the S-wave partin the records and compute the Fourier spectra. Aftersmoothing the spectral amplitudes, the ratios PR4/PR3, PR4!
PR2, and PR4/PR1 are computed. For the linear case, thatis, for the weak motions, we expect to obtain the same ratiosfrom different earthquakes. In reality, the ratios have somevariation basically due to the difference of the incidence angles and the heterogeneity of the medium. The linear transfer
function could be considered as the average of these ratios.
In order to discuss the change of the physical properties dueto the liquefaction, we estimate two transfer functions, forthe records of events before the mainshock and for those ofaftershocks. In addition, the records during the mainshock
are divided into two parts (see the left bottom of Fig. 2a):
one before liquefaction and the other after liquefaction. The
transfer functions during strong motions are estimated fromthe former, and those during liquefaction, from the latter.
Further, the transfer functions during the liquefaction state
are estimated from subsequent aftershocks included in thelong-duration record and one aftershock occurring 3 hr later.
The nonlinear inversion is based on the genetic algorithm (GA) technique in which the evaluation function compares the theoretical and observed transfer functions. In thiscase, the chromosomes are the velocities and damping coefficients for each layer. The performance of the individuals
are evaluated by the normalized residual value obtained from

136 ° 00"

Figure 1. Epicentral location of the earthquakes in Table 2. Biggest star corresponds tothe mainshock.

the fitting of the observed and simulated transfer function by
the revised Haskell-Thompson method with damping factors. After several generations, the performance of the bestindividuals are improved, and the normalized residuals reduced. Then the best individual of all generations with thesmaller residual is taken as the result of the inversion. The
range of velocities for each chromosome is fixed based on
the original velocity structure model obtained from the boring data. The interval between 0 and 83 m is divided into
eight layers with different S-wave velocities according to theoriginal velocity structure. The allowed velocity for eachlayer ranges from 5 to 500 m/sec.
Spectral Analysis
For the analysis in the frequency domain, we select the
S-wave part in the records and compute the Fourier spectra.
After smoothing the spectral amplitudes, the ratios PR4/
PR3, PR4/PR2, and PR4/PR1 are computed. In Figure 3a,
the ratios PR4/PR3 (surface/16 m depth) of the E - W component for the 11 aftershocks have been plotted as an example. The average of those ratios is plotted with the standard deviation in Figure 3b. The linear transfer function
could be considered as the average of the ratios. We estimate
the transfer functions as the averaged spectral ratios between
stations PR4/PR3, PR4/PR2, and PR4/PR1 for the events
recorded before the mainshock and for the aftershocks.
First, we calculate the ratios for E - W and N-S components separately. The averaged ratios for the E - W component are very similar to those for the N-S component for
both before and after the mainshock, as shown in Figure 4.
Therefore, the results of the analysis done for one component
can be generalized for both E - W and N-S components. For
the aftershocks, we consider the 11 events occurring after
the mainshock from 18 January to 14 October 1995. We
exclude the event occurring about 3 hr after the mainshockNonlinearity, Liquefaction, and Velocity Variation of Soft Soil Layers in Port Island, Kobe(a)

Sectior~ A

20

Sectior~ C

10

% 6
o

SeCt~OT~ E
'

2

%

1247

%

-5

Section C

4

2

~

-20

_,oo_

/

,oo
,oo_...........

\/

Mai~shock

6°I

Section B

300~

,oo t

4

-2001 heft i£q.~-300 F f i i i

6

12

(b)

20

24

"<.

28

174 177.seco'ltds

Section A
.

10

20

}_o

o

_,o
-15

°-~°[

-20
,,0

,,,

seconds~

Section H

2O

0
-20

-40_60~ , , , , ,
-80

180

222

I
I
,9~

taa

4

324

second.£

2O

I

°_~o[

"

-4o[ . . . .

318

Section C

Section E%

seconds

284 288 292
seeend.~228

seccq~ds

Section C

5

Sectiony~ 2 o ~=fret aq~--i i i t i

16 20
seco~c]~Section D

\,,o

6O
4O

}°

~o-100~ , o o

,oo

256

264

seco~zds-4°t ....

270\

306

312

1

sec°nds I
300- !...........-':

200 ~

100

}o

' E-:r-:2¥:_2

/

-100
-200

-300
0

i

i

50

100

i /

1 5 0/6

e

i

i

200

250

Ui r--.:!

'::5:2

,\

300

~0

tends
MainshockSectior~ B

300 r
200
I00
0

Sectio~

4

}o

~-100

12

f6

20

seco~cL~

Figure 2.

24

Section F

% 20
"L

0

o

;

-4° f ~ l , , ,-300 F E i h ~ ~ ~ r ~ ,
8

D

28

174 177
secotrcds

180

_-L°

222

Section H

o

-or
228

secend~284 288 292

318

seco'n.ds

Acceleration seismogram (E-W component) recorded during the 1995
Hyogo-ken Nanbu earthquake at (a) surface and (b) 83 m depth.

....

I
324

seco71.ds

1248

J. Aguirre and K. Irikura

(a)

NS
10 2

PR4/PR3

EtF aftershocks

Vs E W (PR4/PRS) BEPORE~'HOCKS............................................
_....,..,.......
,.,.,......,.,.,.....,.

!Z !;!!: ii!i!!i::

10 2

Iii i::i!ii ii:i: i o0li!iiiliii: '. : : : : :

10 ~

NS

Vs EW

(PR4/PR3)

AFTERSHOCKS........

: ! ! ! i: ! !i!: ::: ::: ! ::::::::: .....

:: : :,::::

',

for

::::::::::::::::::::::::::::::::::::::::::::::::i
iii
i ili
:iii!i:;,::!:,2!i~i::::.!:: :::::::::::::::::::::::::': I" :i:iTi-U:,,,:i......

hi!============================================================================================
tO-t

1 0o

101

fO-'
tO°

to'
$

o

10 0

-,i:,.',,"{: '! ~,".. i"_.L.::. _.::... i.. :..!~i. i,.:: ..............
:

10-*

:

:

:

:

i

i

.
:

..............

i-ii !-!:::........ ::.... i_.:...i...L.L::..... :::::::::::::::::::::

i

.

.

;

.

.

.
:

.

:

.

,"IS V s R W ( P R 4 / P R 2 )•~ , - ?

10 ~

:

7
~o t

::

i ........ i...... :::::::::::::::::::::.............

,._• !!!!!!i!!!!!::i!!!!i!!!i:~!:!i!i!?!ii:i!!i~i"':" ":' ':': H ! . . . . . V " : " ":" ":" U:':':: . . . . . . ! ' ":;.i .....................

:N S V~ E ~ g ( P R 4 / p R 2 )t O 2 I .. ..... ..,. . . . . ..... .. .. . . . .. .. ...., .. ..... .. ..,.. ..... .. ..... .. ..... .: . . . . . . . . . . . . . . . . .

................................................
::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::"i'7!'!i!

.....

i'"V

"Vi'!'i'i'~i ......

7":::::::::::::::::::::::::::::::::::::::::::::::::!!i:iS~i!i:?~i;,::!::7:!!}b!!U:!!!}i!!;!!!}ii!

,, .......

-~:i~.~..~.!~!.~!~LZ~:.~.~.~:~.~.~!~:~!~:~!~:~!~.!~!~!~:~:~!~;~....

;....:...;..;..;.

~..1 . . . . . . . . . . . . . .

..... :.,.:.,.:..:..:.:..;
i....i..i..i..Li.-

....

i.

L..!_L.iS..i

:

:

:

i,

i

i i i ii

. . . . . . .

I 0 -2::....;...:...:...:..;: ........

:. . . . . .

:. . . . . . : . . . : . . . :..:..:..:.............

i.......

i. . . . . . i . . . i . .

1.2.

..............

! ........

i......
:

:

i

L,

i

............

....

:. . . . . . . .

,,,i....i.

:

: : :

:

:

.; . . . . . . . . .

; ......

.............

: ......

.............

: .......

: ...........

i..i.i!...i..i..i.i

;

:

: : :

:

:

:

: : :

i

i i iii

.....

10 0

: ......

ito o"~ ,ooo1 ~i:~i!

[ Hz

10 0

PR4/PR8
Z::Z::::::Z::::"'.:i:i'.ZZ:.:ZZZ.,:ZZ:,

EWi:i::":,

.....

tO o

10 t
$

]
,vs

1 o2

tO I
]

l0 2

(b)

.........

tO s

10 ~

frequency''

aftershoclcs

vs

£w

N$

(PR4,/PRI)(pR4/pRt)

Vs EW............................................
::::::::::::::::::::::::::::::::::::::::::::::::::::- ":" ":" i":'::: . . . . . :'" '7 " ' : ' 7 : ":':':: . . . . . : "......,..

.::ZZ: .ZiZ::::Z::Z:.:Z:i:ZZ:i:::Z::::i:Z.ZZZ...:::Z,::Z::.::..............

...

.......

.........

:: :::;:::i::! ::i::!:i: ::::::::: :::::::::: :: :::.:::::: :i .:::::::/:::::U::::::::::::::::::::::! ~ i i ....: i :;:~ : ;: qZil..... i

iiT;il....................... i...... i:ZZ: ::Z : Z:::Z', Z ] I : ' : ; Z!Tiiii...................

Z::ZZZZZ::ZZZ:: Z::Z::: Z;.~Z.%,ZZZZ:Z:ZZ;Zi ::~::::ZZ:Z:~:',iZiZ:ZZ Z Z:Z:XZ:ZZZ:
Io o

to'
$

10 0

..... :....:../,.. ~ .> ~..: ........... ,t::.. :~:..,'v', ~ v';";+:+-"~2~,:'.-:-.'~: ...... ~......

.... :.....i...;:=';.ci.:z-;'.':::::• ~.. ~• .:...:..:• ..:....:..

:

10-I

:

:

:

.......

: ........

: ......

::,,?:F-.i---i.-i..i.............

.. . . . . . . .

:. . . . . .

:. • - i

:

:

: ::

:

i .....

- -:-- -:- --:, -~ •

:

:

101
Y

Figure 4. Spectralratios, from top to bottom PR4/
PR3, PR4/PR2, and PR4/PRI: left, events before themainshock, and right, events after the mainshock. The
continuous line is for the E-W component, and thedashed line is for the N-S component.

: : ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::..... ;...,1...',..;..;;..; . . . . . . . . . . . .
..... ;.. :...:..:..;.:..:..............

:. . . . . . . .
:. . . . . . .

.... :...:...:..:.

ii:.:..:.............

iiiii...........

i .......

:. . . . . . : . - - . ; . . . 1 - . . : - . . : . ; . ; . . . . . . . . . . . . .
; .....
:. . . . . . : . . . . : . . . : . . : . : . . : . : . . . . . . . . . . . . . . . . . . . . .

..... :i...'...ii. i.L! ............. i........ !.......
i-i..i..i.i.l.i: : : : ; :. . . . •. . . .

1 0 - 2",;:,:..':',.:'~ l:,

10 a! ........

: ......

i..

:...:..:..:..:.:

i

i

i

i

.

.............

: .......

............

: ......

i iiii

:

i
lO ~

10 0

frequency[ Hz

I

Figure 3. (a) Spectral ratios PR4/PR3 (E-W component) from 11 aftershocks and (b) the average andstandard deviation for the same events.

because of its discrepancy from the rest of the ratios, probably due to the after effects of the liquefaction as later discussed. The spectral ratios PR4/PR3 for the N-S component
from the events recorded after the mainshock until 19 May
were excluded because of the malfunction of the sensor at
16 m depth (i.e., PR3). Then the transfer function PR4/PR3
for the N - S component from the aftershocks is estimatedusing the records of only two remaining events. Even though
the number of events used for calculation of the transfer
function PR4/PR3 for the N-S component is very small incomparison with the number of events used for the E-W
component, both transfer functions are very similar to eachother.
The averaged spectral ratios accurately represent thetransfer function for the frequency band in which the standard deviation is small. For instance, the standard deviation
for the ratios PR4/PR3 is large at frequencies smaller than
1 Hz and larger than 20 Hz.
Comparing the transfer functions obtained from theevents before the mainshock to the transfer functions fromthe aftershocks in the reliable frequency band (see Fig. 4),
we found some differences; the peak spectral amplitudes before and after the mainshock are located at different frequencies. The change of the averaged ratios is clearer for
PR4fPR3 than that for the rest, PR4/PR2 and PR4/PR1. The

Nonlinearity, Liquefaction, and Velocity Variation of Soft Soil Layers in Port Island, Kobe

peak frequencies for the averaged ratios before the mainshock for the E-W component are 3.2, 3.7, 9.5, and 15.3
Hz, while the peak frequencies for the averaged ratios after
the mainshock are 2.9, 7.4, and 12.8 Hz. All the peaks seem
to shift to lower frequencies (about 20%) for the aftershocks
compared with the peaks for events observed before themainshock. For the same component, the peak frequencies
for the ratios PR4/PR2 are 1.2, 4.5, 6.5, 10.0, and 12.3 Hzfor the events before the mainshock, and 1.2, 4.1, 6.1, and11.3 Hz for the aftershocks. The reduction of the lowest peakfrequency is about 5%, smaller than the PR4/PR3 reduction.
The averaged ratios PR4/PR1 for the same E - W component
show the peaks at 2.0, 4.2, 6.5, 8.0, 11.4, 12.3, and 16.6 Hzfor the events before the mainshock, and 1.9, 3.9, 5.5, 6.3,
9.5, and 11.0 Hz for the aftershocks. The reduction of thelowest peak frequency for this case is also about 5%. Since
these peaks are closely related with the velocity and thickness of the layers, we will try to invert the appropriate parameters that explain the change of peak frequencies for thetransfer functions before and after the mainshock.
Another interesting exploration is done in order to see
the spectral ratio variation in the course of the 360 sec during
the mainshock and subsequent aftershocks (sections A to H
in Fig. 2). We compute the spectral ratios for 11 events including the mainshock and aftershocks successively happening as shown by solid curves in Figure 5. The first two ratios,
M-pI and M-pII, are for the beginning part of the S wavesof the mainshock record before liquefaction, and the part of
S waves after liquefaction. The 3rd to 10th correspond to the
A to H events. The 1 lth is an aftershock about 3 hr after themainshock (exactly 8:58, 17 January 1995). For reference,
the broken curves in Figure 5 show the averaged spectral
ratios for the events before the mainshock and its standard
deviation. The variation of the spectral ratios after the mainshock with respect to the pattern given from the events before the mainshock is dramatic for the ratios PR4/PR3 in thehorizontal components. The records of the A to H events are
mixed by the later phases from the mainshock and the aftershocks just after the mainshock. Nevertheless, we can see
the variations of peak frequencies caused by the rigidity
change of soil layers. For the vertical component, the difference from the reference transfer function is relatively
small for the beginning part of the mainshock before theliquefaction took place and for the 17 January, 8:58 aftershock.
The ratio for the 17 January, 8:58 aftershock for the N S component exhibits an extremely unusual amplification at
low frequencies. The N-S component acceleration record forthis event, recorded at - 16 m itself, has very large amplitudes and looks strange compared with the other component
and the deeper stations. It is very probable that during therecording of this earthquake the sensor enters in the malfunctioned status mentioned before. Therefore, the amplification observed in the N-S component at low frequencies isrelated with an instrumental error.
From the time variation of the spectral ratios in Figure

1249

5, we can see the physical properties changing due to liquefaction and gradually recovering with the lapse time.
Note that this variation after liquefaction is independentof the size of the peak acceleration. The variation in thehorizontal components is larger than that in the vertical component. That is expected since the horizontal components
include mainly S waves that depend on/l (shear modulus),
and the vertical component includes mainly P waves that
depend on a combination of/l and 2 (Lam6 constants). We
consider that the shear moduli of the liquefied layers between the surface and 16 m depth are reduced, close to zero(~ = 0 for liquid), while 2 may not be changed.
Inversion by Genetic Algorithms
Since the transfer function for a given number of horizontal stratified layers is a nonlinear function of the velocity,
thickness, and damping of each layer, the inversion of thevelocities is a nontrivial problem even for the idealized caseof vertical incidence of SH waves. Here, we solve the problem with the help of genetic algorithms, assuming the vertical incidence of mainly SH waves because the horizontal
N-S and E-W components have almost the same tendencies.
The transfer functions are simulated using a revised versionof the Haskell-Thompson method where the damping wasintroduced (Silva, 1976). We assume a constant thickness ofeach layer, according to the original data (Development Bureau, Kobe City, 1991) based on the boring data. The averaged spectral ratios of the events occurring before the mainshock, the mainshock (divided into two parts before and afterliquefaction), the D event in Figure 2, and the aftershocks
are selected to be used to invert the velocity structure by thismethod for the S waves. Only the D event among the A to
H events has amplitude levels large enough compared withlater phases of the mainshock.
For the first generation, several (a given number of individuals per generation) different models are proposed.
Each model or individual has chromosomes chosen from theallowed space randomly. The theoretical transfer functions
computed using the proposed models are compared in turnwith the ratios PR4/PR3, PR4/PR2, and PR4/PR1. The per,formance of each proposed model is judged by the normalized residual. The judgment is used by GA as the parameterto be minimized, and it determines the probability of reproduction, that is, survival, for the next generation. The individuals to produce the next generation are selected based onsuch probability of reproduction. The selected individuals
interchange chromosomes in pairs, producing the new generation. The mutation is considered in order to avoid thetendency to obtain local minimum models. In the mutation,
all the chromosomes have the mutation probability (in ourcase, 0.001) to be replaced by other randomly selected. Obviously, the better individuals will remain through the generations. So, after several generations, the population is improved, and the performance of the better individuals are
better. The goal of producing several generations with sev-

1250

J. Aguirre and K. Irikura

PR4/PR3 NS

10 ~

3

lOt

,oou

PR4/'PR3 E ~

:::

:

:,,.:

: : ::::,ooI
P R 4 / P R 3 UD:,o, ili: ::: i:)~ : : : :.~:-... :i ............

lO-t

1

10o

fO~

10t

!lii! !ii!ii!ii. . . .

~,~,~,.....,,.,

.;,.

0

_' " tI ~*'~ ' 1'~ ~;'*':
0 ' ' 0" ' ~ "

'°'1

100

10~

-.:

:

,-.":"°*,< . . ~ "

lO-t

=1

~.

i i i ii!.iii

:

,

:

M pI

10~

10 0

I0 t ,i: : ::

::

:

M plI

I0 o
10-!

tOo

10~

lO-t

10t

10~

100

! i!

i i;...::i! i iii.i

!

10 o

10° ~

10-t
101

100

lO-t

101

k. . . .

10-t

lO-t

.............

: ~:~ ::::~..;.,., :I.o t o o i ' "~ ,, *~." : : ~'f~':~.". - "~"
lO-t

100

10t

t0 ~

! i:

:

~:~,....~:~.~,o

lO-t

i0-~
tO t

10~

10o

:::

:

~'~],,~ : : :~;.,
:

101

100

:

: ::!:!!!

::::~:~:.:.~:..~

i i ! i!!i!!

i
D

~

'

"
tOt

100

i0 t

.

E

,o-"°°~ ' ~ ": .%#-".':'.:;;..-75~:.~
101

10t

10°

10o
101 . . . . . . . . . . . . . . .

I01

F

101

100
!

101
i

•

tO 0

I 0°

10-I

10-t

;;~i,:...

100
:::

:

:

: :::::;

I0 t

:

101

~

•

100
i!i

: ~.#.fi:;;....'3.~'_-

10 °

10 ~

.o 1oo :"tO-t":.,~ii!iii

i"!~-~,:-':'-'~':.:'.--.~_,%:
i

Io'

~o_,
:1oo
i!i.~. .......
~i i iiiii~
10~

10t

100

!!i

! i i i!!iii

i
H

10°

i'-:~': ~: : :~:i";.":,~, .~, :I0 f/:... ~*~..

10o

too iI :: ~~ :

10~

100

~o'
•

lo'

101

10°

: : ~:~*",:: :::L::: ..... :101 :::

10-t

i

:

10~

100

10~

' i i ~i i iiiii

iii: 5"'~.~::'~L ~ :

10 o

10~

100

10-~

10-t

• ~,,~.'-"~~,~.~t~:.

too ~

;tOt

10-t

t0 ~

~ iiio

:

10t

100
:::

10-t

100

IO-t10° ~100 ~:A

B

fO°~o' .i!!. . . ! ,~:":,~,"
i:.-i!!!i!!
..~-.,~ ~ ~i:..I:,o-,'oo
101:

X,,"; : : : : : :

100

fO ~

• ::

:

:: ; : ' . : :: : ' .

•

10 0

lot

::

: :

10-~
100
10 ~
I0 t . . . . . . . . .

10 ~

tO t

! .;i,~i!i!!ii:!~,'~ . • . . ~-..:~:;

,ooI

100

10~

100

tO~

100
::

~o° ~

10-t

10t
:

~

: :::::::::.........

I0 t
10°

100

10 t

~

Ja,n. 17
8:58

lO-t

¢0o

10t
.r

10-~

;tOo

10~
y

10-t

10o

I0 t
J

Nonlinearity, Liquefaction, and Velocity Variation of Soft Soil Layers in Port Island, Kobe

Figure 5. Variation of spectral ratios for different
sections (in continuous line) with respect to the reference averaged spectral ratios before the mainshockand their standard deviation (dashed lines). Shown atleft, the N-S component, E-W in the center, and UD
in the right. From top to bottom, M-pI (mainshock
before liquefaction), M-pII (mainshock after liquefaction), sections A to H (see Fig. 2), and 17 January,
8:58 aftershock.
eral individuals is to find the best individual that is the solution of the inversion. After several trials, we conclude thatthe population of 24 individuals and 52 generations produces
an appropriate performance of the algorithm. In Figure 6, an
example of the performance of GA for one of the aftershocks
is shown. In this example, the population per generation is24, and the total number of generations is 52, totally, 1248
individuals. One individual corresponds to one simulation.
The genes of each individual for our case are the velocities
and damping coefficients for each layer. To avoid problemsof solutions, physically unreasonable searching space waslimited. The lower velocity limit permitted for each layer is
5 rn/sec, and the upper limit is 500 rrdsec. In this way, therange of velocities is wide enough to allow slower or fastervelocities than the originals, but within a physically reasonable band.
There are several possibilities to solve the problem. We
compromise all the information for the four stations withdifferent modalities. One of them was inverting the velocities and damping factors using the three spectral ratios separately and then interconnecting the results for common
layers. The nonuniqueness of the solution in this kind of1251

problem produces a large dispersion when common layer
characteristics obtained from different ratios were compared.
Other modality consists in minimizing simultaneously thethree spectral ratios PR4/PR3, PR4/PR2, and PR4/PR1. In
this case, the result for each layer was unique, but due to thestrict criterion, the approximation was poor and the residuals
large. Also, we tried to invert interactively velocities with
damping fixed and vice versa. This modality requires an initial model either for damping or for velocities; therefore, theresult is strongly influenced by the initial model selected.
The modality that gives us the best results is to invert the
total structures by parts. We first invert the velocities anddamping factors for 0 to 16 m depth, next for 16 to 32 m
depth, and finally for 32 to 83 m depth using the ratios PR4/
PR3, PR4/PR2, and PR4/PR1, respectively.
For the inversion of the velocities before and after themainshock, the E - W and N-S components are treated independently. First, the velocities and damping factors for thelayers between 0 and 16 m were inverted based on the PR4/
PR3 averaged ratios. Next, using the velocities and damping
factors obtained previously and the average ratios PR4/PR2,
we invert the velocities and damping factors for the interval
16 to 32 m depth. Finally, the velocities between 32 and 83
m depth were inverted using the ratio PR4/PR1, and thevelocities and damping factors obtained in the two previous
inversions for the layers from 0 to 32 m depth. The inverted
velocities for the whole structure are shown in Figures 7 and9. We compare the spectral ratios synthesized using the bestmodels with the corresponding averaged spectral ratios be-

EW component

Vs for---Before

the mainshock

After the mainshock
i

~ x

~

40-

0.6

,

0.4

×x~

I

Depth
(m)

"d

60~

LI

x

x '~1

x

x

x

x
x×

×

x×

~x

Xxx

x

-r

x~ ×

80-

0.2

{

I

]

I

I

I

0

10

20

30
=24

40

50

Oe~e~%l:i,o~,s ( 19e~.

i,~oZi,vi,g~c~Ls)

Figure 6. Example of the performance of genetic
algorithms with 24 individuals per generation and 52
generations.

1 oo,,~,,,,I,,,,L,I,;

¸ ,

ii

~,1',,,,

I,ii,

oo 150 200 250 300 350 400 450 500
S wave velocity (m/s)

Figure 7. S-wave velocity structure inverted fromthe events before the mainshock and from the events
after the mainshock, E-W component.

1252

J. Aguirre and K. Irikura
PR4oPR3 ( 0 - 1 ~ ) E~r bl~J'o~*ePR4oPR8 (0-1~) E W ~Vs for NS componentIO ~!iiiiiii?

!iiil]iii

i!ii!iiii

i]

i iiiiillli!
-Before0. . . . . ?. . ~. . i . m.

I I ....

the mainshock

I ' ' '

After the mainshock
20

10 °

: : '.:::;2:: : :::::::

ltlliHII IIIIIIll

t0-140.1

100

: :
i I
401

40-1

10o

40 f

40 ¸

!

PR4oPR~ (0-83) E ~ before

102 !

Depth
(m)

PR4oPR2 ( 0 - 3 2 ) EW o~ter

10 z

/

60

iiiiiiiii !iiiiilil ill
10 t

I-I
. . . .

w l , ~

i

101

8O
400

10o
; : :::1:::

40-J:','Z:X:::i {IiilI~

;

t0-1

'11

: :

T t iiiiiii
400

I
I
I

i

I

101

10-t

z I

I

10-I

i

i~l

,lll,llll

401

400

.f

i,ii

I,l,,I,,iI,

iI,,i,I....

10ui00'~ 150 200 250 300 3S0 400 450 500y

S wave velocity (m/s)
10.~R4oPR1(0-83) EW Ee,~rsre (rote*ted)

40zPR4oPR1 (0-08) EW ~[ter (rot~tedJ: ::::::::
10 r: ::::::::: ::Figure 9. S-wave velocity structure inverted fromthe events before and after the mainshock, N-S component.

............... ,v. . . . . . . . , . . . . . . . . . . . ...

o

1o o

fO0 . .
; 2222;22

10-t

i

40 -f

i

: ::2222:[HH[[

i ~ LL;m~
100

!

~ I
10 f

,

. . . . . . . . . . . . . . . .

?: iiiiii??
::;:Z::

: 2
10-~

10 -T

~

.

,

,

?: ?iiiiii?:~
i?i: :::::::
'. ::

100

10 t

!

Figure 8. Observed spectral ratios (continuousline) compared with the synthesized spectral ratios(dashed line) using the result of the inversion in Figure 7, for the E-W component.

fore and after the mainshock, for E - W and N - S components
separately in Figures 8 and 10. In Figures 11 and 12, thesimulation of the E - W horizontal motion for an event before
the mainshock and for an aftershock using the corresponding
results of the inversion are compared with the observation,
respectively. In Tables 3 and 4, the results of the inversions
are summarized.
The inversion for the aftershocks produces smaller residuals than the residuals for the events before the mainshock. Using a larger number of layers, this problem may
be reduced, but as the number of parameters increases, thediversity of solutions for the same spectral ratio increases
also.
For interpreting Figures 7 and 9, we should consider
that the velocities inverted from the observed data at nearsurface are better controlled than those at deeper. It is easy
to note since the information about the surficial layers (0 to
16 m depth) is included in all three spectral ratios (PR4/PR3,
PR4/PR2, and PR4/PR1), about the intermediate layers (16to 32 m depth) in two spectral ratios (PR4/PR2 and PR4/
PR1), and about the deeper layers (32 to 83 m depth) in only
one ratio (PR4/PR1). Then, the lack of resolution produces
the dispersion observed from one to another inversion in the
deeper layers, for instance, layers 50 to 61 m depth and 61
to 79 m depth. The velocities inverted for the two layers
located between 19 and 33 m depth from the events before
mainshock for component N - S seems to be affected by thenonuniqueness of solutions. The inversion for the E - W component gives a better result with smaller residuals than thatfor the N - S component. The layer 5 to t6 m depth is the
only one that shows a clear and systematic difference between the models before and after the mainshock in thisanalysis. Since this layer is close to the surface, the parameters are well controlled. Linking these results with otherobservational facts, it is possible to find a physical interpretation, as is discussed later. During the strong ground motions, the S-wave velocity is not a constant but is dependent
on the level of strain, as expected from a nonlinear relationship. But it is possible to obtain an approximation using a
constant velocity considering the concept of "equivalent linear" introduced by Idriss and Seed (1968).
In the case of the mainshock, the inversion is divided
into two parts: pI and pII. The division is to consider apparent velocities before and after the liquefaction occurrence. The liquefaction is assumed to start 15.3 sec after the
beginning of the record (Report of the Japanese Society of
Soil Mechanics, 1996). The S-wave-velocities are found tobe smaller just after the liquefaction occurs and later tend toNonlinearity, Liquefaction, and Velocity Variation of Soft Soil Layers in Port Island, Kobe
PR4oPR,.~ (0-t6) N~ befo~tO =i iiiiiii!
i ii!iiiil!iiiiiii!
i!!!i!ili

:

"

:

"

[ : 7;:::;:::::::::: : ::::::[

:

:

|

Pt74 - O r e....................

lot

10° " ~

Sto~tior~ P1 Do, f,e: 1 9 9 4 / 1 0 / 2 4 11:51 Accelero, f,4.o~. ]~Wl O t'~'rt.ee,'rstrt~'u2~t~r~ (rr~odeL be/ore)

PR4oPR3 (0-t6) NS ~fler

tO •

!i
!!

7
: : ;:::;7',

: ;
6

100

tO-t

8

10

tO-t

10 ~

10 -~

100

/

!

PR4oPR~ (0-32) NS be{ore30 ~i iiiiiiiiiiiiiiilill

10 t

PR4oPR2 (0-32) N,,~ cfter

tO~

12

14

16

sec,

I0 ~

PR3-16rn

]
tO~

1253

~o 1oh~10 t~

Observed^

Sirr~M~ztetl

-10

6

8

Io

12

14

16

sec,

PR2-32"m,

10o

10 o

]

30
20

: : :::::t:

fO-110_t: : ::t:',::

100

: :

fO-t

I0 t

i"i'l l I UlI,,
tO-~

100

$

iI

10 ~
Y

},o
0

-10
Io~R4oPRt (0-83) NS before (rotated)tO t

i iii!!iii

i!ii!iiii!ii: ;::;;::7; :;:;;:S: ;:• - C. -: .[.:.:.: ,3.:...

i
8

Io PR4oPRf (0-83) NS ~.fter (~'otated)

r
10

i
12

i
f4

i
16

sec.
PR 1 - 831"r~t0;. "." ; ; :.:;:....:..:.

1reposed rno~io~oi!ii?'?~tO o•:,iiii:iii
,:.

,:

i i
;,....:..;.

;tO0

~

~

,.

! i!!!!!!!

i ',.'i!!i!i!", ! i

-10
-gO
6

10-I

I0-~

100

10 t

10-t

.f

100

101

8

10

12

14

16

sec.

!

Figure 10. Observed spectral ratios (continuousline) compared with the synthesized spectral ratios(dashed line) using the result of the inversion in Figure 9, for the N-S component.
increase gradually, as was verified with the inversions of thevelocities for the subsequent aftershocks D and the 17 January, 8:58 aftershock occurring 3 hr after the mainshockshown in Figure 13. Even though the level of accelerations
recorded during this time at the surface were small enoughto be considered as a linear response, the spectral ratio analysis shows the nonlinear character of the response. That isthe permanence in the liquefaction state.
After the mainshock, normal ground motions had notbeen recorded until the 18 January 1995, 5:40 aftershock
about 24 hr after the mainshock. After this time, the responseof the layers in the linear range can be predicted by thevelocity structure model obtained previously for the aftershocks.
Nonlinear Simulation
We perform the simulation of the 1995 Hyogo-ken
Nanbu earthquake including the nonlinear characteristics ofFigure 11. Simulation of acceleration waveforms(E-W component) for the event occurring at 11:51 hr
on 24 October 1994 using the velocity structure
model inverted from the events before the mainshock(see Fig. 7).

the soil based on the velocity structure obtained for the
events before the mainshock. This simulation was carriedout by the same program as Perez-Rocha and SanchezSesma (1992), which is a revised version of the Joyner andChen (1975) original program. The nonlinearity is modeledby a Hardin and Druevich (1972) model curve. The convenience of using this model is that the nonlinearity depends
basically on one parameter, the failure strength. The nonlinearity is normalized by the failure strength depending on theproperties of the layer and its position. The layers are represented as a discrete medium in which the stress and displacement are computed carrying out the time integration
step by step using the finite-difference technique. In thisway, the nonlinear stress-strain relation can be introduced
appropriately. Joyner and Chen (1975) presented a detailed
explanation of the method.
The velocities, densities, and failure strength for eachlayer used for the simulation are in Table 5. We use theparameters obtained by the inversion from the events before

1254

J. Aguirre and K. Irikura
Stc~tlo~ 1~1 Do, re: 1 9 9 5 / 0 1 / 2 6 01:01 AcceLe~rzt~.o~ Bgrt B l~nea,r s~muLat~.on (model aftePew)

Table 4
Result of the Inversion of S-Wave Velocities and Damping(Damping shown here)

pR4-Orr~
2O

-to~-

I
8

a

110

~
12

I
14

16

sec.

Obser'ued
20

.

•

~,o
-10

,,,
6

8

1o

I

1

16

sec.

P.R2-..~2rn,
Observed

30 l

l
[

},
6

8

10

12

14

16

sec.
P1~1-83rn!

lm~osed motion

-10 I
-20
6

8

10

12

14

16

Layer

h~ before N-S

hs before E-W

hs after N-S

hs after E-W

135.0
5.1319.0
19.1327.0
27.0-33.0
33.1350.0
50.1361.0
61.1379.0
79.13(83.0)

0.020
0.115
0.057
0.009
0.002
0.105
0.059
0.005

0.011
0.115
0.042
0.011
0.009
0.029
0.012
0.048

0.009
0.115
0.048
0.005
0.011
0.028
0.088
0.012

0.119
0.066
0.078
0.005
0.004
0.003
0.049
0.049

Generally speaking, the simulation fits well the wholeobserved records for the PR2 and PR3 stations. The simulation at the surface fits well the observed signal for the first
15 sec of the observed records in Figures 14 and 15 before
the liquefaction onset. The reason of the misfit after this timeis that our simulations do not consider the water pressure
development, and then it is not able to reproduce the effects
of liquefaction. Nevertheless, it is helpful to estimate thepeak acceleration and to identify the occurrence of liquefaction. Simulations including the water pressure development
have been done by several authors, such as Sato et al. (1996),
Yoshida et al. (1995), Cubrinovski and Ishihara (1996), andElgamal et al. (1996).

Discussion

sec.

Figure 12. Simulation of acceleration waveforms(E-W component) for the event occurring at 01:01 hr
on 26 January 1995 using the velocity structure model
inverted from the events after the mainshock (seeFig. 7).
Table 3Result of the Inversion of S-Wave Velocities and Damping(Velocities shown here)Layer

V s before N-S

Vs before E,-W

Vs after N-S

Vs after E-W

0-5.0
5.1319.0
19.0-27.0
27.0-33.0
33.050.0
50.0--61.0
61.0-79.0
79.13(83.0)

194.0
237.0
191.0
138.0
254.0
208.0
469.0
486.0

159.0
244.0
137.0
169.0
248.0
230.0
443.0
499.0

149.0
195.0
144.0
180.0
244.0
315.0
418.0
478.0

151.0
196.0
136.0
182.0
236.0
235.0
457.0
487.0

the mainshock, and the failure strengths were estimated following Iwasaki (1987). We synthesize the acceleration
waveforms including the nonlinearity at 32, 16, and 0 m
depth imposing the observed motion at 83 m. In Figures 14
and 15, the results of the simulation are compared with theobserved acceleration records.

The horizontal acceleration records in Port Island bythemselves show the influence of the nonlinearity during the
1995 Hyogo-ken Nanbu earthquake since the accelerations
at the surface are reduced in comparison with the accelerations at deeper ground contrary to the amplification expectedfrom linear theory. However, for the UD component, thenonlinearity is not significant before the liquefaction occurrence. Aguirre and Irikura (1995) showed that the observed
acceleration for the UD component is well approximated by
the linear simulation in peak amplitude and general features.
This conclusion is supported by the spectral ratios of Figure
5 in which the nonlinearity of the spectrum during the firstpart of the mainshock (before the liquefaction) for the UDcomponent is smaller than that for the horizontal components.
In the horizontal components, the variation of the averaged spectral ratios after the mainshock in comparison
with the averaged spectral ratios before the mainshock
evinces a change in the S-wave velocity of the layers. Thesechanges related with the peak frequencies are larger for the
averaged ratios PR4/PR3 than for those of PR4/PR2 andPR4/PR1. The change of the S-wave velocities in the surficial layers (0 to 16 m) affects the peak frequencies of theratios PR4/PR2 and PR4/PR1 because they are part of thelayers between surface and - 32 m and between surface and- 83 m. However, the effect of the surficial layers on theNonlinearity, Liquefaction, and Velocity Variation of Soft Soil Layers in Port lsland, Kobe

1255

PR4oPR3 ( 0 - I ~ ) ~gr
10 2

for

Vs

EW

component

J' ......l--'Be°rethe

~

- -i
_ _ I

20

mainshock,, 3-hrs after
the mainshocko: ::::::::: ::2:::::: :;:::;::: ::::::::

: :
: :

i iiiiiiii

i iiiiiiii

i i

i

IOo

lO-f

tO-I

10 °

-I

101

f

40

PR4oPR2 ( 0 - 3 2 ) El#

10z!iiiiiiii i iiiiiiii i i
Depth
(m)
60

i

100

I
I
I

1o-'

: :::::::1: ::::::::: :: ::::::::: ::::::::: ::::::::: : :;:::::

: :
: :

i iiiiiii: :~iiiii:~i:~ i:~
tO-t

10 °

80

1O0

100

200

i i!!iiii? )??)i)).... !,,,,IT,,,I,,~,

300

400

500
r~

S wave velocity (m/s)

lO't

PR4oPR! ( 0 - 8 3 ) EW

102
lll[l~lll~l,lrl,,Ifi iii!!ii i i!
i )?iK? i ii

iO0'1"~i"oi~ii:.

~ :Y:?ill i il!!!!!!!!! !!!!!!i!!' ! ! !
: : ::::=:~ i;;iil;i

lO_t
i0-'l

: ::211;::i ;ii;;iii

10 0

f

Table 5
Parameters Used for the Nonlinear Simulation
Thickness(m)

Density(ton/ms)

S-wave Velocity(m/sec)Damping
FactorFailure Slxength(dyn/crn2)

2.0
3.0
7.6
3.4
3.0
8.0
5.0
1.0
17.0
11.0
18.0
4.0

1.8
1.8
1.8
1.8
1.8
1.5
1.85
1.85
1.85
1.85
1.80
1.90

194.0
194.0
237.0
237.0
237.0
191.0
138.0
138.0
254.0
208.0
469.0
486.0

0.07
0.07
0.08
0.08
0.02
0.50
0.50
0.50
0.50
0.50
0.50
0.50

3.25 E06
5.15 E06
1.21 E07
1.37 E07
1.49 E07
1.64 E07
6.62 E06
6.73 E05
2.84 E06
2.11 E06
11.87 E08
13.86 E08

peak frequencies is much smaller in the ratios PR4/PR2 andPR4/PR1 than in the ratios PR4/PR3. The inverted modelsof S-wave velocity for the 0 to 16 m depth interval obtained
from the average spectral ratio (PR4/PR3) are well distinguished before the mainshock and after the mainshock, because of the clear variation of the peak frequencies in thethose ratios, whereas the ratios PR4/PR1 before and after themainshock do not resolve clearly the velocities for the
deeper layers (from 16 to 83 m depth). However, there is no

: i:
; ~i

10't

Figure 13. S-wave velocity structure inverted from the 8:58 hr, 17 January 1995 eventfor the E-W component in the left, and comparison of observed spectral ratios (continuousline) with the simulated using the inverted velocity structure (dash line) in the right.

indication that the layers below 32 m depth had sufferedliquefaction, and therefore, the possibility of large velocitychanges in deeper layers (>16 m) is smaller than in thesurficial layers ( < 16 m) affected by the liquefaction. In our
analysis, only the second (5 to 16 m) layer presents systematic change with a reduction of 20% of S-wave velocity. The
nonlinear simulations of the acceleration waveforms performed here agrees satisfactorily with the observed accelerations, except for the ones at the surface after the liquefaction occurred. The misfit is natural since our simulations
do not take into account the liquefaction. This resuk meansthat the zone of liquefaction occurrence is restricted to thelayers shallower than 16 m.
Once we fix the second layer (5 to 16 m depth) as thelayer more affected by the liquefaction, let us discuss the
temporal variation of the S-wave velocities for this layer.
Actually, we know that there is a variation of the S-wave
velocity before and after the mainshock caused by the liquefaction. From the linear equivalent velocities inverted
from the pI and plI during the mainshock (before and after
the liquefaction onset), the D aftershock (3 rain after themainshock), and the 8:53 aftershock (3 hr after the mainshock), we can observe a sharp reduction of the S-wave velocity just after the liquefaction onset and gradual recoveryof the S-wave velocity for this second layer. This means that

1256

J. Aguirre and K. Irikura

Station P2~T Date: 1895/01/17 05:46 Accet. E~
Comparison of observations and ID non-linear simulations.
PR4-Om
1500 [% 10oo-obse'rved

Station PRT Date: 1995/01/17 05:46 AceeL NS
Comparison of observatio~zs and ID r~on-linear sirnulatioT~s.
PR4-Om
1500 ~ observe_____~dobserve1000 F~~m~zate~~"~ ~O°t simulated

simulat-500

81500 t
1000}

i

r

i

i,

i

i

h

10

12

14

16
seconds

18

20

22

24

8

10

1500 tobserved

?

i

i

i

i

I

p

10

12

14

18
seco~,,

18

20

obse~e~

.

k

22

24

PR2132rn

i

,

i

i

i

i

i

i

12

14

16
sCqoTzds

18

20

22

24

PR1-83rn

1500

irrqgosed m1000 I

~

q

-500

-500~
i

12

I

14

I

16
seco'n,ds

lt8

210

I

22

l

I

~

i

i

1

i

10

12

14

16
seconds

18

20

22

s

8

0
I

24

~

i

m

i

10

u

t

a

t

e

@

P

i

i

i

i

i

12

14

16
sece'n,eLs

18

20

22

24

Figure 14. Nonlinear simulation for E-W components and the recorded accelerations in Port Island
for the 1995 Hyogo-ken Nanbu earthquake using theparameters in Table 5.

the layer is completely liquefied, and the S-wave velocity is
reduced to be close to zero; that result is in good agreementwith the S-wave velocities obtained for the same interval bySate et al. (1996) and Elgamal (1996). Later, the layer's Swave velocity is increased gradually. We found that 3 hr
after the mainshock the recovery was not complete, but it
was almost completed after 24 hr. Figure 16 shows the process schematically; note that the horizontal axis does nothave a linear scale because the time is elongated and shrankfor different intervals. Solving the differential equation for
1D (one dimension) consolidation model, Terzaghi (1956)
showed that the percentage of consolidation (U) has an exponential dependency with respect to the time { U = 1 exp[-f(t)] }. Considering that the process after liquefaction
involves consolidation, it is expected to follow the same tendency as the curve of consolidation shown by Terzaghi(1956).
Here we discuss the cause of S-wave velocity reductionin the second layer. After the liquefaction, settlement of thesurficial layer (50 cm in average) has been reported on Port

24

~

24

PR1-83m !

1500 t

500~

10

22

observed-500 L

-500~
10

20

18

PR2-32m

1500 t? "°;I

simulated

1000

8% lOOO~1°°° I816
seeond,~

500[

-500

fSO0 t14

observed ~% 1ooo~simulated

8

12

i

8

r

n

p

o

s

e

d

~

t

i

i

i

i

i

i

10

12

14

16

18

20

22

seconds

Figure 15, Nonlinear simulation for N-S component and the recorded accelerations in Port Island forthe 1995 Hyogo-ken Nanbu earthquake using the parameters in Table 5.

300

Before themainshock250

OQOOOI~
After themainshoek

2O0
S wavevelocity 150[ m/s]

Before.~liquefaction ~50tAfter 74 hrs.

- - 0

I I
100

--OOOO000

I
I
•/ ~
] / ~
'

After 3 hrs,

I Ii~-----~After

3 rain,

i'

It

Ti me

Figure 16. Schematic representation of the Swave velocities in Port Island for the second layer (5to 16 m) before, during, and after the 1995 Hyogoken Nanbu earthquake.

24

Nonlinearity, Liquefaction, and Velocity Variation of Soft Soil Layers in Port Island, Kobe
Island (Ishihara et aL, 1996). That indicates that the densityof the second layer increased after the liquefaction. On theother hand, we found in this article a reduction of shear wave
velocity of 20%. Also, the logging data show an S-wave
velocity reduction of about 20% in average for the samelayer (Development Bureau, Kobe City, 1995). A reductionof S-wave velocity in 20% indicates a density increase of
56%, if we consider the/z constant. Considering the density
before the mainshock as 1.8 g/cm 3, the increase of 56% isphysically unreasonable. Then such a change in velocity
only can be explained with a simultaneously increase of density and reduction of/z.
The aftereffects of liquefaction remain at least about 3
hr after the mainshock, as can be seen in the nonlinear response during the 17 January, 08:58 aftershock. The liquefaction causes a permanent variation of the elastic propertiesof the liquefied layers. This implies that the seismic response
at liquefied sites after the earthquake is different from before.

Conclusions
From the spectral analysis, we observed a large variation of the spectral ratios during the strong ground motions
and during the liquefaction state (Fig. 5). The spectral ratios
after the mainshock (and after the liquefaction state) are different from those before the mainshock. The peak frequencies in the spectral ratios after the mainshock are shifted tolower frequencies by about 20%, 5%, and 5% in comparison
to the spectral ratios before the mainshock for the PR4/PR3,
PR4/PR2, and PR4/PR1, respectively.
The results of the S-wave velocity inversions showed
clear differences between the velocity structures obtained
before and after the mainshock; most notably, the S-wave
velocity of the 2nd layer, from 5 to 16 m depth, after themainshock was lower than before.
The biggest change of soil rigidity was observed several
seconds after liquefaction occurrence, then it gradually recovered following an inverse exponential relation with time
(as in the case of consolidation). After the recovery wascompleted, the S-wave velocity of the second layer (5 to 16
m depth) was 20% lower than that before the mainshock.
The other layers did not show a systematic variation.
The aftereffects of liquefaction remain at least 3 hr after
the mainshock but no more than 24 hr.
The nonlinear simulation shows a strong reduction of
the shear modulus during the largest ground shaking before
the liquefaction happened. The observed horizontal peak acceleration was only about 25 % of the peak acceleration expected from the linear theory, while the observed verticalpeak acceleration is very similar to the vertical peak acceleration predicted by the linear theory. The inclusion of thiseffect helps to explain the complex distribution of the peakaccelerations in the Kobe area. It also shows the necessityof considering this effect in the seismic risk analysis.

1257
Acknowledgments

We thank Dr. Tomotaka Iwata for his helpful suggestions. We are also
grateful to L. E. Perez-Roeha for providing the code with which to computethe nonlinear simulation. Detailed comments from David J. Wald and two
anonymous 1~eviewershelped us to improve the manuscript. The acceleration records of the vertical array in Port Island were kindly provided by
Kobe City. The program used for genetic algorithms was GENESIS version
5.0 coded by John J. Grefenstette.

References
Aguirre, J. and K. Irikura (1995). Preliminary analysis of non-linear site
effects at Port Island vertical array station during the 1995 Hyogokennambu earthquake, J. Nat. Disaster Sci. 16, no. 2, 49-58.
Aguirre, J., K. lrikura, and K. Kudo (1994). Estimation of strong groundmotions on hard rock and soft sediment sites in Ashigara valley using
the empirical Green's function method, Bull. Disas. Prey. Res. Inst.,
Kyoto Univ. 44, 45-68.
Ansary, M. A., F. Yamazaki, T. Katayama, and I. Towhata (1995). Analysis
of ground motions at liquefied site during the 1995 Great Hanshin
Earthquake, 23rd Earthquake Eng. Workshop, JSCE, 277-280.
Beresnev, I. A., K.-L. Wen, and Y. T. Yeh (1995). Nonlinear soil amplification: its corroboration in Taiwan, Bull. Seism. Soc. Am. 85, 496-515.
Chin, B.-H. and K. Aki (1991). Simultaneous study of the source, path, andsite effects on strong ground motion during the Loma Prieta earthquake: a preliminary result on pervasive nonlinear site effects, Bull.
Seism. Soc. Am. 81, 1859-1884.
Cubrinovski, M. and K. Ishihara (1996). Assessment of the Kobe Port Island liquefaction through analytical simulation of the vertical arrayrecords, Proc. of the Special Conference on Great Hanshin-Awafi
Earthquake Disasters, Japan Society of Civil Engineers, January, Tokyo, 157-164.
Development Bureau, Kobe City (1991). Report of the soil investigationand installment of strong motion accelerograph (in Japanese).
Development Bureau, Kobe City (1995). Report of the investigation ofground variation in reclaimed land during Hyogo-ken Nanbu earthquake (Port Island, Rokko Island) (in Japanese).
Elgamal, A.-W., M. Zeghal, and E. Parra (1996). Liquefaction of reclaimed
Island in Kobe, Japan, J. Geotech. Eng., ASCE 122, 39-49.
Hardin, B, O. and V. P. Drnevich (1972). Shear rnoduhis and damping insoils: design equations and curves, Proc. Am. Soc. Civil Eng., J. Soil
Mech. Foundations Div. 98, 667-692.
Idriss, I. M. and H. B. Seed (1968). Seismic response of horizontal soil
layers, J. Soil Mech. Foundations Div. ASCE 96.
Ishihara, K., S. Yasuda, and H. Nagase (1996). Soil characteristics andground damage (special issue), Soils Foundations January 1996,
109-118.
Iwasaki, Yoshinori (1987). Dynamic properties of soils. Lecture note oftraining course for the Kansai section of Japanese Geotechhical Society, November 1987, 381-398 (in Japanese).
Joyner, William B. and Albert T. F. Chen (1975). Calculation of nonlinear
ground response in earthquakes, Bull. Seism. Soc. Am. 65, 1315-1336.
Ordaz, M. and E. Faccioli (1994). Site response analysis in the valley of
Mexico: selection of input motion and extent of non-linear soil behaviour, Earthquake Eng. Struct. Dyn. 23, 895-908.
Perez-Rocha, L. E. and F. J. Sanchez-Sesma (1992). Hybrid non-linear
response of soft soils, Proc. lOth World Conf. Earthquake Eng. 2,
973-979.
Report of the Japanese Society of Soil Mechanics (1996). Report on an
investigation of Hanshin-Awaji great earthquake disaster (in Japanese).
Sato, K. T. Kokusho, M. Matsumoto, and E. Yamada (1996). Nonlinear
seismic response and soil property during strong motion (special issue), Soils Foundations January 1996, 41-52.

1258

Satoh, T., T. Sato, and H. Kawase (1995). Nonlinear behavior of soil sediments identified by using borehole records observed at the Ashigara
valley, Japan, Bull. Seism. Soc. Am. 85, 1770-1789.
Silva, W. (1976). Body Waves in a layered anelastic solid, Bull. Seism. Soc.
Am. 66, 1539-1554.
Sugito, M. and H. Kameda (1990). Nonlinear soil amplification model withverification by vertical strong motion array records, Proc. 4th U.S.
Natl. Conf. Earthquake Eng. I, 555-564.
Terzaghi, K. (1956). Theoretical Soil Mechanics, Wiley, New York.
Vucetic, M. and R. Dobry (1991). Effect of soil plasticity on cyclic response, J. Geotech. Eng. ASCE 117, 89-107.
Yoshida, N., S. Nakamura, and I. Suetomi (1995). Nonlinear response of

J. Aguirre and K. Irikura

ground and its evaluation during the 1995 Hyogoken-nanbu earthquake, Proc. of the 23rd Symposium of Earthquake Ground Motion,
Architec. Institute of Japan, Tokyo, 39-52 (in Japanese, with English
abstract).
Disaster Prevention Research InstituteKyoto University
Gokasho, Uji
Kyoto 611, Japan
Manuscript received 13 November 1996.