Journal of Japan Association for Earthquake Engineering, Vol. 10, No. 5, 2010

TEMPORAL VARIATION OF S-WAVE VELOCITY AT
PORT ISLAND VERTICAL ARRAY SITE DURING
AND AFTER THE HYOGOKEN-NANBU
EARTHQUAKE BY USING MODIFIED NORMALIZEDINPUT-OUTPUT MINIMIZATION (MOD-NIOM)METHOD
Santa Man SHRESTHA1, Hideji KAWAKAMI2, Hidenori MOGI3
1

PhD candidate, Dept. of Civil and Environmental Engineering, Saitama University,
Saitama, Japan, santa@kiban.civil.saitama-u.ac.jp
2
Member of JAEE, Professor, Geosphere Research Institute, Saitama University,
Saitama, Japan, kaw@kiban.civil.saitama-u.ac.jp
3
Member of JAEE, Assoc. Prof., Dept. of Civil and Environmental Engineering, Saitama University,
Saitama, Japan, hmogi-2008f@kiban.civil.saitama-u.ac.jp
ABSTRACT: New method of wave propagation analysis is introduced and is used toanalyze the time histories during and after the Hyogoken-nanbu earthquake at Port Island
site in Kobe. The results show large reduction of the S-wave velocity at the surface layer(0-16m) and large reduction of peak amplitude of ground motion at ground surface during
the main shock, indicating liquefaction of that layer. The reflected peaks from groundsurface could not be seen during main shock but can clearly be seen in the aftershocks at
depths. However, the reflected peaks from the boundary of the liquefied and nonliquefied layers can be seen during liquefaction.
Key words: wave propagation, temporal variation of S-wave velocity, liquefaction,
vertical array, amplitude, soil stiffness

INTRODUCTION
Borehole observation (vertical array) records are very useful in the study of local site effects during anearthquake. These records are functions of source characteristics, wave propagation path and the local
site conditions. It is possible to separate the site effect from the source and the path effects using these
vertical array records and can be analyzed to know the behavior of the soil layers. Several methods
and approaches are already in use for wave propagation analysis (Bendat and Piersol, 1971;
Harichandran and Vanmarcke, 1986; Kawakami and Haddadi, 1998; Kramer, 1996; Thomson andDahleh, 1992).
In this paper a new method of wave propagation analysis named as Modified Normalized InputOutput Minimization (Mod-NIOM) is introduced which correlates different time histories observed at- 21 -

various locations on an array. In this method, all the observed time histories are considered as outputs
and they are equally treated. Whereas in the original NIOM method (Kawakami and Haddadi, 1998),
one time history at the ground surface is considered as input and other remaining time histories are
considered as outputs. In NIOM method, only the input is treated specially that causes the input model
similar as a symmetrical single pulse. For this reason, the input model always has to be considered at
the ground surface. A simple pulse signal is considered as an observed input in Mod-NIOM method
for which F    constant. However, if the input model is a rigorous Delta function, the model
response amplitude becomes large at frequencies when calculated spectrum of input motion is almost
zero. In Mod-NIOM method, the sum of the maximum values of the output models is normalized tounity, and considering this as constraint, the sum of the weighted squared Fourier amplitudes isminimized by using Lagrange multiplier method. As all the outputs are dealt equally, the relationshipamong output models obtained from the proposed method is more rational than when they are dealt
unequally without any reason. By using this method, the arrival times and relative peak values forincident and reflected waves can be obtained.
Several studies (Aguirre and Irikura, 1997; Haddadi and Kawakami, 1998; Pavlenko and Irikura,2002) showed that the surface layer of Port Island site was liquefied during the main shock of the
Hyogoken-nanbu earthquake (1995). Aguirre and Irikura (1997) used spectral ratio technique and anonlinear inversion for velocity structure to analyze the data, and their results show the velocity
structure before and after the liquefaction. They observed sharp reduction in S-wave velocity just afterliquefaction and gradual recovery in the layer up to 16m depth. Similarly they observed strong
reduction of the shear modulus during the largest ground shaking before the liquefaction occurred.
Haddadi and Kawakami (1998) used the Normalized Input-Output Minimization (NIOM) method toanalyze the data of the main shock and 5 aftershocks and showed the difference of wave propagation
times during and after liquefaction. In their paper, Haddadi and Kawakami (1998) obtained goodcorrelations among observation records at different seismometer levels and clear peaks in cases ofdifferent aftershocks, but could not find the clear peaks corresponding to wave propagation during themain shock. These studies suggested that surface layer up to 16m depth was liquefied.
In this paper, Mod-NIOM method is used to analyze the wave propagation during and after themain shock of the Hyogoken-nanbu earthquake in the same vertical array site. Correlations amongobservation records at different seismometer levels are obtained even during the main shock by using
Mod-NIOM method. Then the temporal variation of S-wave velocity during and after the main shockis analyzed from the Mod-NIOM results.
METHODOLOGY
Outline of the Mod-NIOM method
Input motion f t  and outputs g l t  for a linear system can be related by

g l t    hl   f t   d

(1a)



and their Fourier transforms can be related byGl i   H l i F i ,

2 

i ,
 i  0,...,N  1, l  1,...,P & i 
N
t 


(1b)

where hl   is the impulse response functions, t the time, and τ the time lag. Similarly, Gl i  are thedifferent outputs of the actual earthquake ground motion records in frequency domain, H l i  thetransfer functions, F i  the input motion in frequency domain,  i the frequency, t the sampling- 22 -

rate in time domain, N the total number of samples, and P the total number of observed time histories
considered as the outputs.
In Mod-NIOM method, all the time histories from the vertical array are considered as „outputs‟. In
Eq. (1a), the actual time histories from the vertical array data are considered as g l t  . A simple pulse
wave is considered as observed „input‟ f t  for which F    constant. The transfer functions H l i 
are always calculated with respect to this input (simple pulse wave) and the actual time histories.
Consider y0 t  and Y0 i  be the input model in time and frequency domain respectively and yl t 

and Yl i  , ( l  1,...,P ) be output models in time and frequency domain respectively of the proposed
Mod-NIOM method. These models in frequency domain can be related similarly as,
Yl i   H l i Y0 i .

(2)

The transfer functions H l i  in Eq. (2) are same as of Eq. (1b), since Y0 i  and Yl i  are the input
and output models of the actual input F i  and outputs Gl i  respectively.
We can write the inverse discrete Fourier transform of the model in general form asymt  

j 2im1 N 1Y i e N ,Nt i 0j   1.

(3)

To formulate the relationship for input-output models, the sum of values of the first output model
y1 (t ) at t m  0 s and the rest of the output models yl (t ) at other times t m  ml t (where l=2,…,P) are1lnormalized to unity which givesP 

y1 0   yl t m  1
l 2

(4a)

l

and transforming Eq. (4a) into frequency domain, we getj 2imP
1 N 1 NY1 i    Yl i e  1.
Nt i 0 
l 2

l

(4b)

Equation (4b) is used as a constraint in the formulation of the input-output models.
To examine the wave propagation in different frequency components (high, medium and low), andto change the spectrum of the models, the derivatives of the Fourier amplitudes are also considered inequation formulation. We use constraint optimization technique to get the optimum models that satisfy
the constraint imposed on them. For this, the sum of weighted squared Fourier amplitudes (amplitudes
of input and output models and their derivatives) is minimized when subject to the constraint of Eq.
(4b). Thus the Lagrange multiplier method gives





N 1
P
2
2
2
2 

2
2
L   c0 Y0  i   k 0 i Y0  i    cl Yl  i   k l  i Yl  i  
i 0 
l 1j 2 im
N

1
P
 1

 
N
 
Y1  i    Yl  i e  1,l 2
  Nt i 0 
l(5)

where L is the Lagrange function,  the Lagrange multiplier, c0, c1,…,cP the weighting constants forcorresponding squared Fourier amplitude, and k0, k1,…,kP the weighting constants for the derivatives- 23 -

of the corresponding Fourier amplitude. The contribution of frequency components is controlled by
the weighting constants k0.
Now substituting Eq. (2) into Eq. (5), we get









N 1
P
2


2
*
2
*
L    c0  k 0 i Y0  i Y0  i    cl  k l  i H l  i  Y0  i Y0  i 
i 0 
l 1

j
2

im
P
 1 N 1 


N
 
 H 1  i    H l  i e
Y0  i   1,l 2
 Nt i 0 

l

(6)

where Y0 i  and Y0* i  are complex conjugate pairs. Since all the output models are dealt withequally, the ratio between ki and ci, (i=0,…,P) are assumed to be constant, then we getN 1 
P

k
2 

2 
*
L   c0   cl H l i  1  0  i Y0 i Y0  i c0
i  0 
l 1




j
2

im
P
 1 N 1 


N
 
 H 1  i    H l  i e
Y0 i   1.
l 2
 Nt i 0 


(7)

l

N
N

N

 1 , i =and   1  i  N  1 ,
2
2
2





Splitting and expanding the terms of Eq. (7) for i = 0, 1  i  converting the expanded terms by using the properties of discrete Fourier transform for real y(t) (i.e.
Y i N   Y i  and Y  N  i   Y *  i  ) and rearranging them gives
N

1

P
P
2 

k 2
k 2
2 
2 


*
*
L   c 0   cl H l  0  1  0  0 Y0  0 Y0  0   2  c 0   cl H l  i  1  0  i Y0  i Y0  i l 1
i 1 l 1
 c 0 c 0




2

P
P
   k  *   2  c 0   cl H l   N  1  0  N Y0   N Y0   N  
H 1  0    H l  0 Y0  0 



c
N

t
l 1
l 2


2 
 2   2  2
0




N
1
2



H l  i e N
  H 1 i   Nt i 1 
l 2
N

P

j 2 im l

1(8)       
Y0  i  
H l   N e jm Y0   N 
 H 1   N   Nt  2  l 2  2    2 Pl j 2imP * 2  **
 H 1  i    H l  i e N Y0  i   .
Nt i 1 
l 2

l

In Eq. (8), Y and Y* are complex conjugates and are of the form Y  x  jy and Y *  x  jy (x and y arereal values), so the value of L may be regarded as a function of either x and y or Y and Y*. It is always
possible to go through the analyses, treating the problem as though x and y are real independent
variables. The minimum condition can be expressed as (Claerbout, 1985)
L
L

 0.
Y Y *

Therefore, minimizing L from Eq. (8) asL 0 yields
Y0 i - 24 -(9)j 2imP
H 1 i    H l i e N
 
*
.
l 2
Y0 i  P
2 Nt  
k0 2  2  c0   cl H l i  1  i  
l 1
 
 c0
 

(10)

 j 2imP
 H 1* i    H l * i e N

 
.
l 2
Y0 i  P
2 Nt  
k0 2  2  c0   cl H l i  1  i  
l 1
 
 c0
 

(11)

l

and

L

Y0 i *

 0 yieldslSimilarly, minimizing L from Eq. (8) as L  0 yields
j 2imP
1 N 1 NH1 i    H l i e
Y0 i   1  0.
Nt i 0 
l 2

l

(12)

Solving Eq. (11) and Eq. (12) in terms of  , we obtain  2Nt 12

 j 2im
P

 *

*
N








H


H

e
H


H

e



 


1
i
l
i
 1 i  l i
N 1 
l 2
l 2






P
i 0
2 
k 2



 c0   cl H l i  1  0 i 
l 1

 c0



j 2im l
N

P

(13)

.

l

Since the values of ml (where l  2,..., P ) are unknown, the optimum value of L is to be foundafter several trials. Substituting the values of Y0* i  , Y0 i  and  from Eqs. (10), (11) and (13)respectively in Eq. (7) and solving, we get
L  Nt 12

 j 2im
P

 *

*
N








H


H

e
H


H

e







1
i
l
i
1
i
l
i
N 1
l 2
l 2






P
i 0
k0 2 2 c0   cl H l  i  1   i 
c0
l 1





P

j 2 im l
N

.

(14)

l

The minimum (optimum) value of L is found out by using Eq. (14) and the corresponding values of
m2 , m3 ,...,mP are fixed regarding optimum value of L. Then the input model can be expressed as
 j 2im
P

 *

*
 Nt  H 1  i    H l  i e N  
l 2

 1
Y0  i   
A
P


k
2

2
0

c0   cl H l  i  1   i  
 
l 1
 c0
 
l

- 25 -(15)

and output models as j 2im
P

 *

*
 Nt H 1  i    H l  i e N  l 2
 1 ,
Yl  i   H l  i 
P
A
k
2 
 c0   cl H l  i  1  0  i 2  


l 1
 c0

 
l

j 2im j 2im
P
P

 *

*
  H 1  i    H l  i e N  H 1  i    H l  i e N  
l 2
l 2

 .
A    Pi 0k0 2 2 c0   cl H l  i  1   i 
l 1

 c0



l

where

(16)

l

N 1

Since we considered all the observed ground motion data as outputs (i.e. Gl i  ), the corresponding
output models are obtained as Yl i  and evaluated by Eq. (16). Therefore, the values of transfer
functions, H l i , (l  1,...,P) , in the input and output models of Eqs. (15) and (16) are calculated from
actual input and output observations by using Eq. (2). These output models in frequency domain canbe converted into simplified output models yl t  in time domain by simply using the inverse Fouriertransforms.
VERIFICATION OF METHOD
To test and verify the effectiveness of the new method, first we apply this method to simulated groundmotion waveforms. To simulate the data (wave forms), a simple time history f0(t) is considered. Other
data f1(t)-f5(t) are simulated as per Eqs. (17a)-(17e) with different amplitudes and different arbitrarily
considered time lags. Some white noises, ni (t ) , (i=1,…,5), are also added. The time histories thus
simulated as
f1 (t )  0.1 f 0 (t )  n1 (t )f 2 (t )  f 0 (t  0.30)  f 0 (t )  n2 (t )f 3 (t )  f 0 (t  0.45)  f 0 (t  0.15)  n3 (t )f 4 (t )  f 0 (t  0.65)  f 0 (t  0.35)  n4 (t )f 5 (t )  f 0 (t  0.75)  f 0 (t  0.45)  n5 (t )

(17a)
(17b)
(17c)
(17d)
(17e)

These simulated time histories are shown in Fig. 1(a). These are analyzed with the proposed
method taking weighting coefficients of c0  1, k0  0.001 and cl  1, (l  1) . These values are taken to
obtain good spectrum of the models and also to analyze these time histories in medium frequency
components. The results (output models y1(t)-y5(t)) are shown in Fig. 1(b). The simple peaks on the
result correspond to the arrival times (shifting times) of the simulated waves. The difference of arrival
times between data levels is the same as the shifting time used in Eqs. (17a)-(17e). Similarly the
relative peak amplitudes obtained are similar to the relative amplitude of the actual simulated time
histories. In Fig. 1(b), the reflected peaks (peaks in positive time domain) appear smaller than theincident peaks (peaks in negative time domain). Since the sum of values of the incident clear peaks ofall the output models is normalized to unity in this method, the effect of different noise components
added in each simulated time histories made the reflected peaks appear smaller than the incident peaks.

- 26 -

(a)(b)
Fig. 1 (a) Simulated time histories as per Eqs. (17a)-(17e), (b) Mod-NIOM results for simulated data
Comparing from the peak at surface or peak of 1st output model, the peaks appearing in the left
hand side correspond to the incident waves, whereas the peaks appearing in the right hand sidecorrespond to the reflected waves. As the wave traverse from depth to the surface, the peakcorresponding to incident wave appears initially in the output model at lower level as seen in Fig. 1(b).
Similarly, the peak corresponding to reflected wave can be seen (if any) later in the model at lower
level. The difference of time of the clear peaks in the output models at two levels gives the one way
travel time of the wave between these two levels. In this way, the propagation times of the wave are
obtained from the Mod-NIOM analysis. Once the propagation times between two levels are known,
the S-wave velocity in that layer (between two levels) can be found as
v

H
,
T

(18)

Here H is the height (thickness) of the layer and T the S-wave propagation time from bottom to top
level of this layer.
ANALYSIS OF ACTUAL TIME HISTORIES AT PORT ISLAND VERTICAL ARRAY SITE
Port Island is a manmade island near the Kobe city. It suffered heavily and was liquefied during themain shock of the Hyogoken-nanbu earthquake also known as Kobe earthquake on 17th January, 1995.
The vertical array site is located in the northwest corner of Port Island. The Port Island vertical arraysite consists of four seismometers placed vertically at 0m, 16m, 32m and 83m depth from the groundsurface. Various earthquake events during and after the main shock of the Hyogoken-nanbu
earthquake are analyzed by using Mod-NIOM method. The list of earthquake events used in theanalysis is summarized in Table 1. To find the wave propagation pattern and variation of S-wave
velocities during the main shock, the horizontal components of the main shock are analyzed with amoving time window of 5.12s from 10s to 70s in 2.5s intervals. The first rectangle (a) in Fig. 2 shows
the total analyzed time from 10s to 70s. Similarly, at later parts (aftershocks) as shown by rectangles(b)-(e) in later times in Fig. 2, four time windows of 40s data are analyzed from 160s to 360s. The
horizontal components of 0-70s data of the main shock, at different levels of the array, are shown in
Fig.3. The soil profile and other relevant soil properties at Port Island site is given in Table 2(Department of Urban Development, Kobe city, 1995). Although the soil profile consists of several
layers of soil between two seismometers, in this analysis, only one layer is assumed between two
seismometers. Therefore, three layers namely 1st layer (surface layer i.e. from 0-16m depth), 2nd layer(intermediate layer i.e. 16-32m depth) and 3rd layer (base layer i.e. 32-83m depth) are considered for- 27 -

the analysis. The reference S-wave velocities for these three layers are calculated by taking an averageof S-wave velocities from the down-hole boring data. After averaging from the down-hole data, the Swave velocities in surface, intermediate and base layer are taken as 196m/s, 203m/s and 313m/s
respectively.

Fig. 2 Time history of EW component of main shock at 83m depth of Port Island vertical array site
during the Hyogoken-nanbu earthquake (Department of Urban Development, Kobe city, 1995)
Box (a) shows the total analyzed time from 10s to 70s during which moving windows of
5.12s at an interval of 2.5s are considered, Boxes (b)-(e) represent time windows of 40s data

Fig. 3 Time histories of horizontal components (0-70s) of main shock observed at Port Island verticalarray site during the Hyogoken-nanbu earthquake of 17th Jan, 1995 (Department of Urban
Development, Kobe city, 1995)- 28 -

Table 1 List of earthquake events used in this analysis(Department of Urban Development, Kobe city, 1995)
Date

Time

1995.01.17

5:46

1995.01.17

5:53

1995.01.17

Epicenter

Magnitude

LatitudeLongitude

Depth (km)N34 ﾟ 35.7'

E135 ﾟ 02.2'

16

7.2

8:58

N34 ﾟ 34.9'

E135 ﾟ 00.5'

19

4.5

1995.01.18

5:25

N34 ﾟ 41.5'

E135 ﾟ 11.1'

15

4.3

1995.01.19

5:10

N34 ﾟ 40.8'

E135 ﾟ 09.7'

10

3.2

1995.01.26

1:01

N34 ﾟ 46.0'

E135 ﾟ 16.0'

12

3.3

1995.01.28

8:06

1995.01.29

16:02

N34 ﾟ 41.0'

E135 ﾟ 10.8'

14

3.6

1995.02.02

16:19

N34 ﾟ 41.4'

E135 ﾟ 08.5'

18

4.1

1995.02.08

15:40

1995.02.18

21:37

N34 ﾟ 26.2'

E134 ﾟ 48.9'

16

4.8

1995.02.24

8:03

N34 ﾟ 42.7'

E135 ﾟ 12.4'

14

3.3

1995.03.05

10:04

N34 ﾟ 44.2'

E135 ﾟ 14.4'

13

3.2

1995.04.06

10:50

N34 ﾟ 47.4'

E135 ﾟ 19.2'

12

4

1995.05.04

5:53

N34 ﾟ 41.7'

E135 ﾟ 11.1'

15

3.6

1995.05.08

2:36

N34 ﾟ 42.6'

E135 ﾟ 12.8'

14

3.3

1995.05.15

7:33

N34 ﾟ 38.8'

E135 ﾟ 07.9'

13

3.4

1995.05.28

10:34

N34 ﾟ 40.4'

E135 ﾟ 10.1'

14

3.2

1995.06.19

8:38

N34 ﾟ 45.8'

E135 ﾟ 16.2'

13

3.4

1995.07.24

4:00

N34 ﾟ 55.4'

E135 ﾟ 33.1'

14

4.1

1995.09.12

6:30

N34 ﾟ 41.9'

E135 ﾟ 11.9'

13

3.9

1995.10.14

2:04

N34 ﾟ 37.6'

E135 ﾟ 06.8'

15

4.5

1995.11.13

15:31

1995.12.27

9:18

N34 ﾟ 26.5'

E134 ﾟ 47.7'

10

3.8

1996.01.25

20:17

N34 ﾟ 46.4'

E135 ﾟ 20.7'

15

3.3

Table 2 Soil profile and other relevant soil properties at Port Island site (Department of Urban
Development, Kobe city, 1995)

Depth (m)Soil type0 - 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
50.0 - 61.0
61.0 - 79.0
79.0 - 85.0

Sandy gravel
Sandy gravel
Reclaimed
Sandy gravel
Sand with gravel
Alluvial clay
Diluvial sand
Sand and clay
Sand with gravelinter-layered
stratum
Diluvial sandDiluvial clay
Sand with gravel- 29 -

P-wavevelocity(km/s)
0.26
0.33
0.78
1.48
1.18
1.33
1.53
1.61
1.61
2.00

S-wavevelocity(km/s)
0.17
0.21
0.21
0.18
0.25
0.31
0.35
0.30
0.32

Poisson
Ratio
0.13
0.32
0.46
0.49
0.49
0.48
0.48
0.48
0.48
0.49

Location ofseismometer(m)GL-0.0GL-16.0
GL-32.0

GL-83.0

RESULTS AND DISCUSSIONS
Fig. 4 shows some of the analysis results for different time windows (moving time windows of 5.12s
from 10s to 70s in an interval of 2.5s) of EW components of the main shock of the Hyogoken-nanbu
earthquake (Fig. 4a-4c) and EW components of different aftershocks (Fig. 4d and 4e) as indicated in
the top left side of each plot. Similarly, Fig. 5 shows the results of NS components. In these results(Figs. 4a-4c and 5a-5c), the peak amplitude of output model at the surface is very small in all the casesof the main shock due to reduced soil stiffness and increased damping during liquefaction. But thepeak at the surface is amplified in case of all the aftershocks as seen in Figs. 4d-4e and 5d-5e,
indicating the recovery from liquefaction.
Fig. 6 shows the results for different time windows of the EW components of the main shock
records and the aftershocks. Similarly, Fig. 7 shows the results of NS components. Figs. 6a and 7ashow the results of the first 20s duration of the main shock, Figs. 6b and 7b show the results of all time
windows of the main shock after the liquefaction, and Figs. 6c and 7c show the results of all theaftershocks. In Figs. 6a-6b and 7a-7b, the peaks at 32m depth are referenced to 0s and thecorresponding peaks at other levels are plotted with respect to this reference. Due to this reason, thepeaks in the 0m and 16m layers are seen in the positive side of the time. In case of Figs. 6c and 7c, thepeaks at ground surface are referenced to 0s for all the aftershocks analyzed. These referencing aredone to show the propagation of reflected waves more clearly during and after the liquefaction.
From these figures (Figs. 4-7), the reflected peaks from the surface could not be seen in the 16m
depth but the reflected peaks from the boundary of liquefied and non-liquefied layers can be seen atthe 32m depth during the liquefaction (as marked by red arrows in Figs. 6b and 7b). This means thetime histories at 32m depth consist of reflected waves from the boundary of liquefied and nonliquefied layers during the liquefaction in addition to incident waves. From the results, the boundary of
the liquefied and non-liquefied layers supposed to lie near or some distance above 16m depth. Graphs
plotted for aftershocks clearly show the reflected peaks from the surface at lower depths (as marked by
red arrows in Figs. 6c and 7c). This also indicates the recovery from liquefaction. In Figs 6c and 7c,
the reflected peaks on positive time domain are smaller than incident peaks. This is not the effect ofliquefaction. In case of liquefaction, we could not find any peak (waves) reflected from the ground in
the liquefied layer (Figs. 6a-6b and 7a-7b). The decrease in the reflected peaks compared to incident
peaks are due to the noise associated in the individual time histories analyzed. The similar effect wasobserved by Kawakami and Haddadi (1998) in their paper for a strong motion records at a TTRL and
Chiba sites. These sites didn‟t show nonlinearity during the analyzed strong motions. However, theresults from NIOM analysis showed decreased reflected peaks compared to incident peaks.
Fig. 8 shows the temporal change of S-wave velocities during and after the main shock of the
Hyogoken-nanbu earthquake. Since the data is not open to public, and we didn‟t have access to thetime histories before the main shock, the records before the main shock could not be analyzed.
Therefore for comparison purpose, we chose the weighted average velocities calculated from thedown-hole boring data. These weighted average velocities for three layers are shown by dotted lines.
Solid rectangles (red), circles (blue) and triangles (pink) represent the S-wave velocities of EW
components of the main shock and aftershocks at the surface, intermediate and base layer respectively,
obtained from Mod-NIOM analysis. Similarly, open symbols represent the S-wave velocities of NS
components in these three layers. Fig. 8a shows the S-wave velocities for different time windows ofhorizontal components of the main shock. In Fig. 8a, the horizontal axis represents the elapsed times
representing the middles (centers) of the time windows used in the analysis. In this figure, the S-wave
velocity in the surface layer is significantly reduced during 17-40s of the main shock. Then it
gradually increases with time. Similarly, the S-wave velocity in the intermediate layer is also reduced
during 12-15s, but increases gradually after that. The S-wave velocity in the base layer is slightly
reduced during first 30s of the main shock and then it increases afterwards.

- 30 -

yi(t)

a)

yi(t)

b)

yi(t)

c)

yi(t)

d)

yi(t)

e)

Fig. 4 Results for EW components of main shock (different time windows) and aftershocks. The time
windows for main shock and the date of aftershocks are indicated in the top left corner of eachgraph.

- 31 -

yi(t)

a)

yi(t)

b)

yi(t)

c)

yi(t)

d)

yi(t)

e)

Fig. 5 Results for NS components of main shock (different time windows) and aftershocks. The time
windows for main shock and the date of aftershocks are indicated in the top left corner of eachgraph.

- 32 -

yi(t)

(a)

yi(t)

(b)

yi(t)

(c)

Fig. 6 Analysis results plotted together for EW components of main shock (top and middle) andaftershocks (bottom). For main shock, graphs are plotted from 12.5-20.0s (top) and from 22.570.0s time windows (middle). Red arrows indicate the reflected peaks.

- 33 -

yi(t)

(a)

yi(t)

(b)

yi(t)

(c)

Fig. 7 Analysis results plotted together for NS components of main shock (top and middle) andaftershocks (bottom). For main shock, graphs are plotted from 12.5-20.0s (top) and from 22.570.0s time windows (middle). Red arrows indicate the reflected peaks.

- 34 -

(a)

(b)

(c)
Fig. 8 Temporal change of S-wave velocities during a) the main shock, b) the aftershocks
up to 2 days and c) rest of the aftershocks. The dotted lines represent the
averaged velocities calculated from the down hole measurements for three layers.
At around 40-70s (of the EW component) and on the later parts of the aftershocks (after 200 days), theS-wave velocities at the base seem to have some discontinuities. This is because of the large depth orlarge velocity in the base layer. For single time increment, the difference of S-wave velocity in thebase layer is around 21m/s where as in the surface layer; it is around 2m/s.
Fig. 8b shows the S-wave velocities for aftershock up to 2days (48hrs). The horizontal axis is theelapsed time in hours. Fig. 8c shows the S-wave velocities for other aftershocks up to 1 year after themain shock. The horizontal axis is the elapsed time after the main shock in days. In Fig. 8c, it can beseen that the S-wave velocities in base layer remain almost same for all the aftershocks whereas theincrease in velocity in surface layer seems to be continuing at a slower pace. The temporal change of
the S-wave velocity is the effect of both a) nonlinearity due to the strong ground motion and b) therecovery of the material property as a function of elapsed time. The sharp reduction of S-wave
velocity is due to the liquefaction which reduces the stiffness and increase the damping. The material- 35 -

properties gradually get recovered on passage of time after the liquefaction which increases the Swave velocity. The temporal change of the S-wave velocity includes the change of the stiffness andthe change of the damping. However, the more dominant, of these two reasons, is the change ofstiffness of the surface layer during the main shock.
From Fig. 8c, it can be seen that the analyzed S-wave velocities for base layer is larger than the
averaged velocity calculated from down-hole data. It also shows that the S-wave velocities calculated
by Mod-NIOM method are consistent for all the horizontal components of the aftershocks. Since theanalyzed velocities in this layer are consistent throughout all the aftershocks, the averaged velocity
calculated from down-hole data for this layer seem to have some errors. The errors in the down-holedata may be due to the effects of material and radiation damping on waveforms which may cause
difficulty in identification of the S-wave arrivals at depths greater than 30-60m during the down-holetests (Krammer, 1996). The other reason of errors in down-hole data may be due to the errors inreading the starting times of the time histories at different receivers during the analysis by using leastsquare method.
From analysis, the velocity at the end of liquefied state at surface layer (0-16m) is about 30% of theaveraged velocity obtained from down-hole data for that layer. Although full recovery of S-wave
velocity didn‟t take place after the end of liquefied state, the recovery process seems to have continued
for about a year.
CONCLUSIONS
A new simplified method of wave propagation analysis named as Mod-NIOM is developed which iscapable of showing the correlation between the observed time histories
The analysis results of the horizontal components of the main shock of the Hyogoken-nanbu
earthquake at Port Island site show large reduction in the amplitude of peak at ground surface
indicating liquefaction. Due to this, the reflected peaks in downward direction from ground surface
could not be seen but the reflected peaks from the boundary of the liquefied and non-liquefied layers
can be seen during the main shock (liquefied duration).
The analysis results of the aftershocks show increased peak amplitudes at the ground surface
indicating regaining of soil stiffness and reduction in damping. Due to this, the reflected peaks fromground surface are clearly visible at lower depths.
Temporal variations of S-wave velocities show that the S-wave velocity in surface layer decreased
sharply during the principal motion of the main shock and it gradually increased after that. The
recovery in velocity in this layer still continued for a year later. The S-wave velocity in the base layerslightly reduced during the main shock but recovered to a consistent value in the aftershocks.
ACKNOWLEDGEMENT
The authors would like to thank the Department of Urban Development, Kobe City for providing theearthquake observation records at Port Island vertical array site.
REFERENCES
Aguirre J. and Irikura K. (1997). “Nonlinearity, liquefaction, and velocity variation of soft soil layers
in Port Island, Kobe, during the Hyogo-ken Nanbu earthquake,” Bulletin of the SeismologicalSociety of America, Vol. 87, No. 5, 1244-1258.
Bendat J. S. and Piersol A. G. (1971). Random Data: Analysis and Measurement Procedures, John
Wiley & Sons.
Claerbout J. F. (1985). Fundamentals of Geophysical Data Processing, Blackwell Scientific
Publications, 109.

- 36 -

Department of Urban Development, Kobe city (1995). Report on ground deformation of reclaimed
land due to Kobe earthquake (in Japanese).
Haddadi H.R. and Kawakami H. (1998). “Effect of liquefaction on ground motion during the
Hyogoken-nanbu earthquake, 1995, in Japan by using NIOM method,” The Effects of Surface
Geology on Seismic Motion, Irikura, Kudo, Okada & Sasatani (eds), Balkema, Rotterdam, Vol. 2,1015-1022.
Harichandran R. S. and Vanmarcke E. H. (1986). “Stochastic variation of earthquake ground motionsin space and time,” Journal of Engineering Mechanics, Vol. 112, No. 2, ASCE, 154 -174.
Kawakami H. and Haddadi H. R. (1998). “Modeling wave propagation by using Normalized InputOutput Minimization (NIOM),” Soil Dynamics and Earthquake Engineering, Vol. 17, 117-126.
Kramer S. L. (1996). Geotechnical Earthquake Engineering, Prentice Hall, New Jersey.
Pavlenko O. and Irikura K. (2002). “Changes in shear moduli of liquefied and nonliquefied soils
during the 1995 Kobe earthquake and its aftershocks at three vertical-array sites,” Bulletin of the
Seismological Society of America, Vol. 92, No. 5, 1952-1969.
Thomson W. T. and Dahleh M. D. (1992). Theory of Vibration with Applications, Prentice Hall, New
Jersey.

(Submitted: February 12, 2010)(Accepted: September 01, 2010)- 37 -