3D waveform simlation in Kobe of the 1995 Hyogoken-Nanbu earthquake by FDM using with discontinuousgrids
S. AoiNational Research Institute for Earth Science and Disaster Prevention
H. Sekiguchi, T. Iwata
Disaster Prevention Research Institute, Kyoto University
H. Fujiwara
National Research Institute for Earth Science and Disaster Prevention

ABSTRACT : We performed a simulation of the 1995 Hyogoken-Nanbu Earthquake with 3D
finite-difference method (FDM) using discontinuous grids. In the waveform simulation using FDM, so longas we use the grids with uniform spacing for the calculation, the grid spacing is determined by the shortest
wave length in the 3D simulation model space. In this study, we used discontinuous grids that consist ofsmall grids (spacing of 50 m) in the region shallower than 2.6 km where there are low-velocity sedimentarylayers, and coarse grids in the deeper region. In this way, we were able to reduce both the time and memory
required for the computation to approximately 1/5 of what was needed in a calculation using small, uniform
grids in the entire region.

1 INTRODUCTION
In Kobe, the most heavy damage of the 1995
Hyogoken-Nanbu Earthquake was distributed on along zone and this is called the “Damaged Belt”.
This particular distribution is due to theunderground structure and the propagation ofruptures. Underground surveys using a reflection
method or a refraction method were carried out
after the earthquake, and the relatively detailed
underground structure has been revealed (e.g.
Huzita 1996). We also learned the details of thesource process through inversions using data suchas strong motion waveforms and geodetic data (e.g.
Yoshida et al. 1996, Sekigughi et al. 1998). Basedon such information, 3D simulations wereperformed (e.g. Furumura & Koketsu 1998, Iwataet al. 1998, Kawase et al. 1998, Pitarka et al. 1998).
These simulations underestimate the amplitude,
although they do successfully explain thedistribution of maximum amplitude qualitatively.
One of the reasons is that their insufficientconsideration of the low-velocity superficial layer(e.g. Iwata et al. 1998).
As long as we simulate the waveforms with
the FDM using uniform grids, the grid spacing isdetermined by the shortest wave length in themodel space. Hence, though the low-velocity
layer exists only near the surface, we are still
S. Aoi, H. Sekiguchi, T. Iwata & H. Fujiwara-1-obliged to use small grids. Thus it was difficult to
introduce the low-velocity layer into the model,
because of the drastic increasing of the grid number,in other words, required memory and computation
time. Aoi & Fujiwara (1998a, b) formulated FDM
using discontinuous grids that consist of small grids
only near the surface where the low-velocity
sedimentary layers exist, and coarser grids for thedeeper region where the wave velocity is higher.
This has shown to be very effective in reducing
computing requirements. In this study, we applied
this method to the 1995 Hyogoken-Nanbu
Earthquake. A superficial layer with the S-wave
velocity of 400 m/s, which is realistic according to
the structure surveys, was introduced into the modelin this study.

2 METHOD
We use discontinuous grids (Aoi & Fujiwara
1998b) that consist of staggered grids with two
different grid spacing (Fig. 1). In Region I, thegrid spacing is small (spacing h1 is 50 m in thiscase), whereas in Region II, the spacing is threetimes coarser (spacing is h2=3h1). These tworegions overlap in depth by 3h1 /2. As to the termsof partial derivatives inside Regions I and II, we use
the fourth-order finite-difference approximation

Top of the Region IIAABCDE：Vx：Vy：Vz：τii：τxy：τxz：τyz

ERegion IIRegion IEABCD

Bottom of the Region IFig.1 : (center) 3D discontinuous grid system. (left) Two transections on the top and at the bottomof the overlapping region of Regions I and II, where the elimination or the insertion of grids arenecessary. (right) Two profiles of the discontinuous grids.
condition of Clayton and Engquist (1977).

(e.g. Levander 1988), denoted as
1f i = { c 0 ( f i +1 / 2 − f i −1 / 2 ) − c1 ( f i + 3 / 2 − f i −3 / 2 )}
h(1).

5

S. Aoi, H. Sekiguchi, T. Iwata & H. Fujiwara-2-8Kmntio

areulaimes34 40’

for the term of partial derivatives.
In theconnecting region between the two regions, we cancalculate neither the velocities nor the stresses
using
thesecond-orderfinite-difference
approximation (the left side of Fig. 1). Hence,these values are calculated by the linear
interpolation and the elimination from the values ofthe other region. We should pay attention to thefact that we perform the interpolation and theelimination of all the variables within horizontal
planes where the grid points concerned exist.
In order to minimize the artificial reflected
waves from the boundary of calculation region, weuse both the non-reflecting boundary condition of
Cerjan et al. (1985) and the absorbing boundary

aActive fault(2),rctu

3Du
str

20 km

135 00’34 50’

Intensity VII(JMA-scale)Fault model we usedkm

1f i = ( f i +1 / 2 − f i −1 / 2 )
h.6
13

Here
wecannotuse
thefourth-orderfinite-difference approximation for the term ofpartial derivatives concerning z in the adjacentregion of the connecting region between Regions Iand II. In such a region, we use the second-orderfinite-difference approximation (e.g. Virieux 1984),

3 UNDERGROUND STRUCTURE MODEL

34 30’（where c0 = 9 / 8 、 c1 = 1 / 24 ）135 10’135 20’135 30’Fig. 2 : Map of Kobe area. The rectangle
shows the area we are modeling in the study.
The location of the heavily damaged zone isalso shown [intensity VII, JMA scale, digitizedby Koketsu (1997)]. The solid line shows thefault segment locations, and the dotted line
shows active faults, traced by Huzita &Kazama 1983.

02010Table 2. Locations of stations.

50403050

04030
201010Fig. 3 : Bedrock topography in the Kobe area used inthe 3D modeling. The shadow line shows the fault
segment locations.

Staion
code
KBU
MOT
SKU
PIS
RKI
FKA
JMA
KB1
KB2
KB3
KB4

Latitude(N)
34.72500
34.72500
34.73056
34.67000
34.68877
34.69500
34.68833
34.72034.71234.70434.698Longitude(E)
135.24000
135.28100
135.25000
135.20800
135.27248
135.21083
135.17944
135.247135.254135.260
135.265

Station name
Kobe Univ.
Motoyama
Shinkobe Substation
Port Island Borehole
Rokko Island
FukiaiJMA Kobe Station
Stations for SS projectSKH
KBUKB1 MOT
KB2
KB3
KB4
FKARKI34 40’JMA34 45’

10 km

PIS
Intensity VII(JMA-scale)Fault segments we used
Observation point
Stations for SS project

135 10’135 15’135 20’Fig. 4 : Map of observation points and thepoints designated by this simultaneous
simulation project. The locations of fault
segments and the heavily damaged zone arealso shown.
Table 1. Model parameters.

Rock layer

Sedimentary basin

Model A
Thickness = 50 %Vp = 2.0 km/s
Vs = 0.6 km/s
Dens. = 1.8 g/cm3
Q = 80Model B
Thickness = 16.7 %Vp = 1.8 km/s, Vs = 0.4 km/s
Dens. = 1.75 g/cm3, Q = 60
Thickness = 33.3 %Vp = 2.2 km/s, Vs = 0.8 km/s
Dens. = 1.85 g/cm3, Q = 100

Thickness = 50 % of sedimentary layer
Vp = 2.4 km/s
Vs = 1.1 km/s
Q = 120
Density = 2.1 g/cm3
Depth = 〜 4 km depth
Vp = 5.5 km/s
Vs = 3.2 km/s
Density = 2.6 g/cm3
Q = 300
Depth = 4 km 〜 17.8 km depth
Vp = 6.0 km/s
Vs = 3.46 km/s
Q = 400
Density = 2.7 g/cm3
Depth = 17.8 km depth 〜
Vp = 6.7 km/s
Vs = 3.87 km/s
Q = 500
Density = 2.8 g/cm3(1998) modified the bedrock topography of Pitarka
et al. (1998), and divided its sedimentary layer into
two layers of the same depth. We call the modelof Iwata et al. (1998) Model A. The origin and thestrike of the model are (34.584N, 134.881E) andN58E, respectively, and its size is 58.0 km (N58E),
13.6 Km (N148E), and 23.0 km (depth) (Fig. 2).
Fig. 3 shows the bedrock topography of our models.
The velocity structure of rock layer shown in Table
1 is from a structure model used for the hypocenter
determination in this area by the RCEP-DPRI,
Kyoto University. According to Iwata et al.
(1998), although the use of the Model A and thesource model of Sekiguchi et al. (1998) led to a
successful
representation
ofwaveformcharacteristics, the amplitudes are underestimated.
They point out that this is due to the effects ofamplification caused by the superficial layer with
low S-wave velocity, which is not taken into
account sufficiently.
In this study, we use a model B where theupper sedimentary layer of the two layers of Model
A is divided into two by 1:2 (Table 1). The
S-wave velocity in the first and the second layers ofModel B are 400 m/s and 800 m/s, respectively.
This means that the interval velocity of these layers
is identical to the S-wave velocity of the first layerof Model A, 600 m/s.

4 RESULT
4.1 Comparison with the observed waveforms

In Kobe area, the underground structure has been
investigated in detail using reflection and refraction
methods and boring experiments (e.g. Huzita 1996).
Using a bedrock shape estimated from thisinformation, Pitarka et al. (1998) performed asimulation of 3D FDM using a model withhomogeneous sedimentary layer. Iwata et al.
S. Aoi, H. Sekiguchi, T. Iwata & H. Fujiwara-3-Fig. 4 shows the observation points where the dataof the mainshock are available, the points
designated by this simultaneous simulation project,
and the areas severely damaged where the JMA
intensity was VII (digitized by Koketsu 1997).
KOB, KBU and SKH are the bedrock sites and therest of observation points are located on theFN (N148E)[m/s]
1.0FP (N58E)

UDObs.
Model A

0.0

Model B

-1.0

No obs. data

0.0

15.0 [s]

30.0

0.0

15.0 [s]

30.0

0.0

15.0 [s]30.0Fig. 5 : Comparisons between observed waveforms and synthetics of fault normal (N148E), fault
parallel(N58E) and up-down components by Models A and B (0.1-1.0 Hz). The numbers on theright side of each trace shows the maximum amplitude in m/s.
sedimentary layer (Table 2). Fig. 5 shows theobserved waveforms (0.1-1.0 Hz) of theseobservation points and the synthesized waveformsof Models A and B. The horizontal component ofRKI and the fault parallel component of FKA wereunderestimated on the stations on the sedimentarylayer in the simulation using Model A. In thesimulation using Model B, the amplitude increased
by approximately 30% and 70%, respectively,
reproducing better the maximum amplitude of theobserved waveforms. The bedrock sites such as
SKH and JMA that gave the underestimated results
logically do not show any improvement. On thesestations, the weathering of the very superficial
S. Aoi, H. Sekiguchi, T. Iwata & H. Fujiwara-4-granite might be causing the reduction of S-wave
velocity. It is difficult to introduce such a thinlayer into models even with the grids of 50 m that
we used. The dominant periods of PR1 are notclose to the observed waveform. This might bedue to the effects of the very thin and softsuperficial layer, and this is also difficult tointroduce into models for the same reason.
4.2 Synthetics on KBU-RKI lineFig. 6 shows the velocity waveforms (0.1Hz-1Hz)of two horizontal components as well as themaximum amplitude on each station of Models A

FP
RKI
KB4

KB2
KB1

Model A
Model B

Sediment

KB3

Model ARockKBU

Model B

FN
RKI
KB4

KB2
Model A
Model B

KB1

5.0

10.0

20.0

T [s]

30.0

Model B
5.0

10.0

20.0

T [s]

30.0

0.0

Model ARockKBU

1.5
1.0
[m/s]
0.5

Sediment

KB3

Fig. 6 : Synthetics of fault parallel and fault normal components by Models A and B on KBU-RKIline (0.1 - 1.0 Hz). The right side shows the maximum amplitude on each station by both models.
and B on KBU-RKI line. The waveforms of thetwo models on the rock sites match very closely,
while the stations on the sediment show that themaximum amplitude of Model B is greater than thatof Model A by approximately 20 – 100%.
KBU is on the boundary of the rock and thesediment, and the depth of the basement reaches
1km within 1km in south-east. In the zone
between 1 km and 3 km from the rock, themaximum amplitude of fault parallel component byModel B exceeds 0.8 m/s. In addition to theamplification effect of soft sedimentary layer, this isalso due to basin-edge effects. Slightly less than 1
km further in the south-east, the fault normalcomponent exceeds 0.8 m/s, and still closer to thecoast, the amplitude decreases gradually. The
S. Aoi, H. Sekiguchi, T. Iwata & H. Fujiwara-5-

width of the zone where the maximum amplitude offault normal component exceeds 0.8 m/s isapproximately 2.5 km.
This zone of large
amplitude corresponds to the “Damaged Belt”.

5 DISCUSSIONS AND CONCLUSIONS
We simulated strong ground motions in the sourcearea during the 1995 Hyogoken-Nanbu Earthquakeby fourth-order 3D FDM using discontinuous grids.
Compared to the computation for same period
range (0.1-1.0 Hz) with uniform grids, thediscontinuous grids reduced the computational
requirements to approximately 1/5. This leads tosmall grid spacing of 50 m in FDM calculation forsedimentary layer.
Using a structure model
consisting of three sedimentary layers with theS-wave velocity of 400 m/s, 800 m/s and 1100 m/s,
and a source model (Sekiguchi et al. 1998) which isapplicable up to 1Hz, we simulated strong groundmotions in the source area and compared with theobserved waveforms on several stations. The
maximum amplitudes of the horizontal components
of RKI by Model B are larger by approximately
30% than those by Model A, they become almost
the same level of those of the observed waveforms.
At the stations on the sedimentary layer on theKBU-RKI line, we found that the maximum
amplitudes by Model B against those by Model A
are approximately 20-100% larger.
A thin superficial layer of less than 100-400 m
thickness with the low S-wave velocity of 400 m/s
influences even waveforms with relatively long
period up to approximately 1Hz. The influence is
not limited to the amplification by the low-velocity
layer: it also influences the secondary generationand the propagation of the surface waves.
Therefore, the influence of the superficial layer
cannot be evaluated using 1D response functions.
The existence of such superficial layers isconfirmed by structure surveys, and their influence
has to be taken into account when we synthesize
waveforms even of relatively long period. The
FDM using discontinuous grids is a very useful toolfor a simulation considering the superficial layer oflow velocity.
ACKNOWLEDGEMENT
Computation time was provided by the Supercomputer centerof National Research Institute for Earth Science and Disaster
Prevention. The data used here are collected and distributedby the Japanese Working Group on Effects of Surface
Geology (ESG) on Seismic Motion, Association for
Earthquake Disaster Prevention, for the Kobe Simultaneous
Simulation Project during the second InternationalSymposium on ESG 98. We used GMT Version 3 forfigures.

REFERENCES
Aoi, S. & H. Fujiwara 1998a. 3-D finite-difference
method using discontinuous grids , Bull. Seism.
Soc. Am. (in printing).
Aoi, S. & H. Fujiwara 1998b. 3-D Fourth-Order
Finite-Difference
Seismograms
Using
DisontinuousGrids,
Japan
Earthquake
Engineering Symposium Proceedings (inS. Aoi, H. Sekiguchi, T. Iwata & H. Fujiwara-6-printing).
Cerjan, C., D. Kosloff, R. Kosloff, & M. Reshef
1985. A nonreflecting boundary condition for
discrete acoustic & elastic wave equations,Geophysics 50: 705-708.
Clayton, R. & B. Engquist 1977. Absorbing
boundary conditions for acoustic & elastic waveequations. Bull. Seism. Soc. Am. 67: 1529-1540.
Furumura, T. & K. Koketsu 1998.
Specific
distribution of ground motion during the 1995
Kobe earthquake and its generation mechanism,
Geophysical Research Letters, 25: 785-788.
Huzita, K. 1996. On survey results of active faults
in Osaka-Kobe area, Proc. 9th seminar on
studying active faults, ‘on deep structure of
Osaka-bay area’, 1-1 (in Japanese).
Huzita, K. & T. Kazama 1983. Geology of theKobe district. Quadrangle Series Scale 1:50,000,
Geological Survey of Japan, Tsukuba, 115p (inJapanese).
Iwata, T., H. Sekiguchi, A. Pitarka, K. Kamae & K.
Irikura 1998. Evaluation of Strong Ground
Motions in the Source Area, druing the 1995
Hyogoken-Nanbu (Kobe) earthquake, Japan
Earthquake Engineering Symposium Proceedings(in printing).
Kawase, H., S. Matsushima, R.W. Graves & P.G.
Somerville 1998. Strong Motion Simulation
inKobe with Special Reference to the Damage
Belt Formation during the Hyogo-ken Nanbu
Earthquake of 1995, Japan Earthquake
Engineering Symposium Proceedings (inprinting).
Koketsu. K. 1998. An overview of the rupture
process, Report on the Hanshin-Awaji
Earthquake Disaster, Common Series, 2: 168-75.
Tokyo: Maruzen (in Japanese).
Levandar,
A.
R.
1988.
Fourth-orderfinite-difference P-SV seismograms, Geophysics
53: 1425-1436.
Pitarka, A., K. Irikura, T. Iwata & H. Sekiguchi
1998.
Three-dimensional simulation of thenear-fault ground motion for the 1995
Hyogoken-nanbu, Japan, earthquake, Bull. Seism.
Soc. Am. 88, 428-440.
Virieux, J. 1986. P-SV wave propagation inheterogeneous
media:
Velocity-stressfinite-difference method, Geophysics 51:
889-901.
Yoshida, S., K. Koketsu, B. Shibazaki, T. Sagiya &Y. Yoshida 1996. Joint inversion of near- andfar-field waveforms and geodetic data for therupture process of the 1995 Kobe earthquake, J.

Phys. Earth, 44, 437-454

S. Aoi, H. Sekiguchi, T. Iwata & H. Fujiwara-7-