Bulletinofthe SeismologicalSocietyofAmerma,Vol.76, No. 6,pp. 1623-1640,December1986

R U P T U R E PROCESS OF THE 1983 JAPAN SEA (AKITA-OKI)
EARTHQUAKE USING A WAVEFORM INVERSION METHODBY EIICHI FUKUYAMA AND KOJIRO IRIKURAABSTRACTThe rupture process of the 1983 Japan Sea earthquake (MAMA ---- 7.7) isdetermined by applying waveform inversion of displacement-type, strong motion
records based on the Bayesian method. The synthesis of the seismic motions
from the main shock is made using the records from the two aftershocks (MAMA= 7.1 and MaMA = 6.1) as empirical Green's functions. The main shock fault
consists of two subfaults. The aftershock (M = 6.1) occurring south of the main
shock fault is used for the southern subfault empirical Green's function and theaftershock (M = 7.1) occurring north of the main shock fault is used for thenorthern subfault. In this synthesis, it is not necessary to calculate propagation
path and local site effects. Therefore, the information about the source process
can be extracted from the comparison between observed seismic motions andthe synthetic ones without detailed information on the velocity structure along
the whole path from the source to each station. The inversion for the rupture
process is attempted for both a line source model and an area source model. In
both models, the main shock fault surface is divided into several elements whose
sizes are determined from the scaling relations between the main shock and its
aftershocks. The slip displacement and the rupture starting time on each element
are estimated as model parameters in the inversion. A good convergence isobtained after about 10 iteration steps for the line source case and after 5
iteration steps for the area source case, The results are summarized: (1) the
average rupture velocity is approximately 2.5 km/sec for the southern subfault
and 2.0 km/sec for the northern subfault, and (2) there is a large slip displacement
near the rupture initiation point and the north edge of the fault. Thus, there isconsiderable heterogeneity in rupture propagation and slip distribution over thefault plane.

INTRODUCTION
High quality records of strong ground motions and detailed studies of synthetic
seismograms have demonstrated that fault rupture during large earthquakes isnonuniform both in space and time. The heterogeneity of the fault rupture process
has important consequences not only for the nature of the earthquake source but
also for the prediction of strong ground motions during large earthquakes. The
purpose of this paper is to show how the irregularity of fault motions can be derived
by inverting near-field seismic records.
Heaton and Helmberger (1979), who determined the slip distribution over thefault surface of the 1971 San Fernando earthquake, demonstrated the heterogeneityof the fault slip by comparing synthetic seismograms based on the Cagniard-de
Hoop solution with the strong ground motions observed during the earthquake by
the forward modeling. Since then, several authors have applied inversion techniques
to estimate the fault heterogeneity or fault motion complexity. The techniques
employed in these studies include three different types. The first one consists ofdetermining the dislocation at each grid point of a finely segmented net that covers
the fault, by the point-by-point inversion method (e.g., Hartzell and Heaton, 1983),
The second consists of calculating moment rate functions at several stations
1623

1624

EIICHI FUKUYAMA AND KOJIRO IRIKURA

surrounding the source area using a deconvolution method to determine where thelarge moment was released on the fault plane (e.g., Kikuchi and Kanamori, 1982;
Ruff and Kanamori, 1983). The third consists of determining the nonuniformnessof dislocation and rupture propagation over the fault directly by using the isochron
method (e.g., Spudich and Frazer, 1984).
To examine the source process, we need to separate the source effect from theothers, such as local site and path effects. All of these techniques discussed above
require some knowledge of the geological structure along the propagation path toeach recording station. Recently, a method for synthesizing seismograms from larger
events has been successfully developed using the seismograms of smaller events as
empirical Green's functions (Hartzell, 1978; Kanamori, 1979; Hadley and Helmberger, 1980; Irikura and Muramatsu, 1982). This method has the advantage that weneed not remove the propagation effects. Then, we propose this '~burth method to
extract the inhomogeneity of the source process. First, we calculate synthetic
seismograms for a main shock using aftershock records as empirical Green's functions. Second, we calculate the cross-correlation error functions between observed
and synthesized seismograms. Third, by applying the Bayesian method for anonlinear inversion problem, we obtain the rupture propagation and slip distribution
over the fault, which minimize the cross-correlation error functions. Some of theadvantages of this method are: (1) we need not take into account the path effects;(2) we need not know the absolute time in each seismogram; and (3) we can directly
combine a p r i o r i information as initial values to the inversion procedure. The
features of (1) and (2) increase the amount of seismic data available for invertingfor the source process.
We applied this method to the 1983 Japan Sea (Akita-Oki) earthquake. This
event occurred on 26 May 1983 in the Japan Sea off the western coast of northern
Honshu [Japan Meteorological Agency (JMA), 1984]. The JMA network recorded
strong motions with low-gain displacement seismographs of both the main shockand the two largest aftershocks at stations located at distances between 150 and
300 km away from the epicenter surrounding the source region in a semi-circle. In
this paper, we have applied the new method discussed above to quantify theirregularities of the rupture process on the fault by inverting these data. We also
discuss the resolution and reliability of the above method.
SYNTHETIC SEISMOGRAMS USING EMPIRICAL GREEN'S FUNCTIONS

To calculate synthetic seismograms using aftershock records as empirical Green's
functions, we follow the method by Irikura (1983). First, we assume the similarity
relation between large and small earthquakes (Aki, 1967). From the similarity
relation, we can derive the formula below (Kanamori and Anderson, 1975; Geller,1976).
L / L e = W / W e = D / D e = t/te = ( M o / M o e ) 1/3(1)

where L and W are the length and width of the fault, D is the final offset of the slip
displacement, t is the rise time, Mo is the seismic moment, and the parameters with
subscript "e" are for the small event used as an empirical Green's function. (Those
without subscript "e" are for the large events to be synthesized.) We divide the fault
surface into several elements. When assuming a line source model, we divide thefault surface into N elements, and when assuming an area source model, we divide

RUPTURE PROCESS OF THE 1983 JAPAN SEA EARTHQUAKE

1625

the fault surface into N x N elements where N is determined by
N = (Mo/Mo~) 1/S.

(2)

We calculate synthetic seismograms by the equation (3) derived by Irikura (1983),
which gives the ground motion at the observation station as
U ( x , t) =

N

N

~

~

R,/Re.(rJr,) p m , . N . Ue(x, t - t~)(3)~=1 h ~ lfor the line source model, andN 2NU(x, t) =Z R,/Re. (re~r,) p m,. U , ( x , t - t~)
z=l k=lfor the area source model, where(4)t,~ = (5 - re)lYe + t,h

and R, is the radiation pattern, r, is the focal distance for the ith element, and Vc isthe phase velocity of the wave of interest, p is 1 for the body waves and ½ for thesurface waves. The usage of subscript "e" in equations (3) and (4) is for the small
event as in equation (1). Here, we assume thatt,k = t , +

(k-

1)t~

(5)

where ti is the starting time of the rupture at ith element, and te is the rise time ofthe small event. A schematic illustration for the synthesis is shown in Figure 1.
We introduce a new parameter m, which is not defined in Irikura (1983). If thescaling relation (1) holds exactly at all elements of the fault, m, should be unity.
However, in reality, the scaling relation is considered to give the average information
about the fault, and so there must be some local stress release fluctuations around
Observation

Station~ !

i TM e l e m e n t

Aftershock- •Mainshock

FIG. 1. Schematic model which shows the relation between the aftershock and the ith element of themain shock.

1626

EIICHI FUKUYAMA AND KOJ1RO IRIKURA

the average. Now, for the sake of simplification, we here assume the rise time isconstant at each element. Therefore, the fluctuations affect only the parameter ms,i.e., m~ stands for the inhomogeneity over the fault. We estimate parameters m, andt~ in the inversion procedure to be introduced in the next section.

INVERSION PROCEDURE

The parameters to be estimated by inverting observed data are m, and ti inequation (3). Although parameter m, can be linearly linked to the synthetic seismograms, parameter t, cannot be, and then we need to use a nonlinear inversion
technique. We use here the cross-correlation error function proposed by Mellman(1980) as the data used in the inversion procedure. This approach has the advantagethat we can use the information from the seismic time series even if we do notaccurately know the arrival times of the seismic waves. The cross-correlation error
function y~ for the ith station isy~ = 1 -f

O , ( t ) S d t + tm~x) d t~(IfO,(t) 2 dtJ[.js~(t) 2 d t(6)where O,(t) is the observed seismogram for the ith station, and S i ( t ) is thesynthesized seismogram, tm~x is the lag time which minimizes yi. According to theBayesian approach (Jackson and Matsu'ura, 1985), we can easily find the solution
by the implicit equation
A TE-le + D-ld = 0(7)where A, e, d, E, and D are defined in equations (All), (A4), (A5), (A6), and (A7)in the Appendix, respectively. We outline the derivation of equation (7) in the
Appendix (see also Jackson and Matsu'ura, 1985).
We solve the nonlinear problem (7) using the Gauss-Newton algorithm derived
by Tarantola and Valette (1982). The (k + 1)th solution xk+l can be obtained by
xk+l = xk + b M k - l r k(8)

M k = A k T E - 1 A k + D -1

(9)

rk = A k T E - l e k + D - l d k(lO)

where b is a constant between 0 and 1. We also consider the partitioned resolution
matrices to separate the estimators operating on a p r i o r i information and observational data from each other (Jackson, 1979). These resolution matrices ( H A for
observational data and K I for a p r i o r i data) are defined as
H = M - l A T E -1
K = M - 1 D -1
M = ATE-1A+ D -1(11)

and
I = HA + KI

(12)

1627

RUPTURE PROCESS OF THE 1983 JAPAN SEA EARTHQUAKE

where I is the unit matrix. (HA)tr and (KI)tr are the traces (sum of the diagonal
elements) of HA and KI, respectively, and the ratio (HA)tr/(KI)tr indicates how
much information we have derived from the observed data.
The algorithm we use here to find the solution of equation (7) consists of two
steps. Step 1: we estimate the parameters through equations (8), (9), and (10). The
convergence criterion of step 1 is that rh of equation (10) is either equal or very
close to zero. Step 2: we calculate (HA)tr/(KI)t~ and compare this ratio with theconvergence criterion (we adopt 0.1 for the present analysis for the line source
model) of step 2. If the ratio is greater than the criterion, we change the initial
parameters for the final parameters of step 1, and then we perform step 1 algorithm
repeatedly until (HA)tr/(KI)tr is less than the criterion.
DATA

We synthesize the ground motions from the main shock of the 1983 Japan Sea
earthquake using the records from the two largest aftershocks (M = 7.1 and M =6.1) and compare the synthetic seismograms with the observed seismograms fromthe main shock. The locations of the JMA network stations we use in this analysis
and the aftershock distribution are shown in Figure 2. The epicenters of the main
shock and the two aftershocks are shown as solid circles, and their focal mechanisms
are also shown to the left. The epicentral distances to these stations range from 150
to 300 km, and the stations are distributed encircling the source area with thecoverage angle of approximately 190 °.
We use the two horizontal component records (NS and EW) with JMA-type
displacement strong motion seismographs. At each station, two types of the strong
motion seismographs with different magnification are installed, i.e., one has unity
magnification (S-type) and the other, 100 times magnification (D-type). The response characteristics of these seismographs are shown in Figure 3. We use therecords with unity magnification for the main shock and the largest aftershock at
all stations and for the second largest aftershock at HAC, SEN, and NII stations
and the records with 100 times magnification for the second largest aftershock at
SUT, HAK, OFU, and AIK stations. Since these data were recorded on paper, wefirst digitized the records, corrected for the instrument response, interpolated at

L2o', I, 50 krn ',LHAK
MainshockI ~Fault - - ~

41=

/
AIK~

L~HAC
~
OFU
SEN

~JNll
40°

200 km
I

!38.

13g,

I

140-

FIG. 2. (Left) The locations of the main shock and two largest aftershocks we used in this analysis
are shown by the solid circles. Focal mechanisms of these events are shown to the left (after Dziewonski
et al., 1983) with lower hemisphere, equal-area projection and "+" denotes the pressure axis. (Right) The
locations of JMA stations we used in this analysis.

1628

E I I C H I F U K U Y A M A AND K O J I R O I R I K U R ARESPONSE

100

L

z
0

i

CURVE
i

P

10

100

D

S

0,1

I
0.1

I
PERIOD(sec)

FIG. 3. Characteristic response curves of the JMA strong motion seismometers with magnificatioJof 1 and natural period of 6 sec (S-type), and seismometers with magnification of 100 and a natura
period of 5 sec (D-type). Both seismometers have the same damping factor (h = 0.55).

13g40

FAULT
s,XIS
~XIS
)f 1983Japan
Earthquake

139 00139 200

20 krnFIG. 4. Geological structure of the main shock region (after Maritime Safety Agency, 1984).

0.1-sec intervals, and then bandpass-filtered from 0.05 to 1.0 Hz. The length of eaci
seismic record for the present analysis is 1 min, which contains the major portion~,with large amplitudes.
FAULT MODEL
Initially, we considered a single fault plane model, and we calculated syntheti(
seismograms using only one aftershock (M = 6.1). This oversimplified mode:
did not yield an acceptable correspondence between synthetics and observec
seismograms.
Next, we used a double fault plane model consisting of two fault surfaces wit[
different strikes as proposed already by several authors, e.g., Shimazaki and Mor:
(1983) from surface wave analysis. Geologic mapping (Figure 4) also suggests that
there are two different fault groups with different strike directions in this region
As the rupture process would likely be different on each fault surface (Sato, 1983
1985; Ishikawa et al., 1984), we used the aftershock (M = 7.1) and the aftershock(M = 6.1) for the northern subfault and the southern subfault as empirical Green'~functions, respectively. Source parameters and respective fault dimensions adopted
here are shown in Table 1. We divided each fault surface into N subfaults (linesource case) or N × N subfaults (area source case), where N is estimated from

RUPTURE PROCESS OF THE 1983 JAPAN SEA EARTHQUAKE

1629

equation (2). T h e seismic m o m e n t ratio between the m a i n shock a n d the corresponding aftershock is e s t i m a t e d by calculating the spectral ratios in the low-frequency
range at stations S U T a n d S E N . T h e ratios are
M o ( M = 7 . 7 ) / M o ( M = 6.1) - 300

andM o ( M = 7 . 7 ) / M o ( M = 7.1) - 10.

Following Mori a n d S h i m a z a k i (1983), we assume t h a t two-thirds of the m o m e n t is
distributed on the n o r t h e r n fault surface a n d one third on the southern fault surface.
T a k i n g this into account for each subfault, N in equation (2) is calculated as follows
N (southern subfault)= [Mo(southern subfault of m a i n s h o c k ) / M o ( M = 6.1)] 1/3 = 5(13)

N (northern subfault)= [Mo(northern subfault of m a i n s h o c k ) / M o ( M = 7.1)] 1/3 = 2.

(14)

F r o m the above estimates, the m a i n shock fault surface is segmented as shown in
Figure 5.
T h e aftershock area dimensions of the two events used as empirical G r e e n ' sfunctions (M = 6.1 a n d M = 7.1) are a p p r o x i m a t e l y 6 x 6 and 20 x 20 kin,
TABLE 1
SOURCE PARAMETERS FOR THE FAULT M O D E L

LengthWidth
StrikeDip(a)NSouthern Subfault

60 km
30 km
N15°E
20°

30 km
30 km
N10°W
20°
N

(b)

10°

~

Northern Subfault

10° ~"°#2 rn
Okm

30k

,1/

t6

2 0 °~ 1 2Line Source Model

Area Source Model

FIG. 5. Fault model assumed in this analysis: (a) line source model and (b) area source model.

1630

EIICH[ FUKUYAMA AND KOJIRO IRIKURA

respectively (JMA, 1984). These sizes are coincident with the element dimensiom
we determined from the moment ratio. This suggests the scaling relation on th(fault length or width to be approximately held among these three events (M = 7.7
7.1, and 6.1) used in this analysis.
Although parameters m~ and t, are considered to be model parameters, the ris(time at each subfault is assumed to be constant. The observed seismograms havepredominant frequency less than 0.2 Hz. Then, we can assume that the deviatiorof the rise time has a smaller influence on the main motions compared with theother parameters. In calculating synthetic seismograms, we assume the rise time o:the main shock is 4.0 sec so that the rise times of the aftershocks (M = 6.1 and= 7.1) are estimated to be 0.8 and 2.0 sec, respectively, from equations (1), (2), (13)and (14).
Several authors have reported that the second subevent occurred 23 to 26 se(after the first subevent occurrence (Mori and Shimazaki, 1983; Sato, 1983, 1985
Ishikawa e t al., 1984). We adopt 25 sec for the retarded time of the second subevent
occurrence when constructing the initial model.
Sato (1983, 1985) showed that the rupture propagated parallel to the fault strikefrom south to north from the analysis of accelerograms obtained at 10 statiom,
located near the fault region (50 to 200 km). Therefore, we first attempt to apply
finite line source model to estimate the source process as shown in Figure 5a. But
at the station perpendicular to the strike direction, the waveforms could not be fully
explained using the line source model. Next, we tried to estimate the source processusing a finite area source model (Figure 5b).
RESULTS OF THE INVERSION

L i n e s o u r c e case. We need a starting model, i.e., a set of initial values to solve
the nonlinear inversion problem. We assume that the horizontal component seismograms used in this analysis consist mainly of S wave, and we perform a forward
test to get approximate estimates of the parameters. Initially, we assumed thatparameters rn~ are all unity and t~ are given as the rupture propagated uniformly
with the velocity Vr. We varied the rupture velocity Vr to find parameters tz whichmakes synthetics agree reasonably well with observations. From the results of theanalysis of teleseismic data (Mori and Skimazaki, 1983; Ishikawa e t al., 1984) andnear-field accelerograms (Sato, 1983, 1985), we assumed that the rupture initiatedat element 2, shown in Figure 5a, propagated bilaterally in the southern subfault(elements 1 to 5) and restarted at element 6 in the northern subfault 25 sec after
the initial starting time. The attempt to determine the initial values was made b)
forward modeling, and the resulting synthetic seismograms are shown in Figure 6
This result suggests a rupture velocity of 2.0 km/sec. Thus, we choose parameter t
for each subfault consistent with the rupture propagation of 2.0 km/sec as thestarting model (Table 2). We assume the covariance matrices E and D defined b)(A6) and (A7) are as follows: all off-diagonal elements of both E and D are equal t(
zero. All diagonal elements of E are equal to 0.3. Diagonal elements of D corresponding to the relative moment releases of northern and southern subfaults and rupture
starting times of northern and southern subfaults are equal to 2.0, 1.0, 4.0, and 2.0
respectively.
After approximately ten iteration steps, we obtained a final solution that satisfied
the convergence criterion [ ( H A ) t r / ( K I ) t ~ = 0.1] of the inversion. The starting mode]and the final model are shown in Figure 7 and Table 2. Standard deviations fol

RUPTURE PROCESS OF THE 1983 JAPAN SEA EARTHQUAKEFORWARD

1631

TEST

su T E[

W

S

Vr = 1.5 km/sec

0 729
Vr :

1.75

0.807
V r = 2.0

0.920
Vr : 2.25
0.886

~

.

~

~

"

V r = 2.5

0.697
V- = 2,75
0.495

0.476

,J

15sec
.

.

-.

.

.

.

Observed
Synthesized.

FI6. 6. Forward test for the line source model. Rupture velocity V, is the only variable. Solid lines
are observed seismograms, and dotted lines are synthesized. The numbers plotted below the waves are
the values of the cross-correlation.
TABLE 2
RESULTS FOR THE LINE SOURCE MODEL
ElementNo
1

2

3

4

5

6

7

Starting model
t,*
S.D.t
m,*
S.D.

3.0
1.4
1.0
1.0

0.0
1.4
1.0
1.0

3.0
1.4
1.0
1.0

6.0
1.4
1.0
1.0

9.0
1.4
1.0
1.0

25.0
2.0
1.0
1.4

40.0
2.0
1.0
1.4

Final model
t,
S.D.
m~
S.D.

2.5
0.55
1.3
0.17

0.5
0.39
1.0
0.16

3.1
0.56
1.0
0.14

5.4
0.54
0.8
0.13

7.4
0.50
0.6
0.14

24.8
0.64
0.7
0.09

40.4
0.60
1.3
0.26

* t, and m, represent the rupture starting time and relative moment
release, respectively.
¢ S.D. represents the standard deviation.

model parameters are also shown in Table 2. The waveforms which are calculated
from the final fault model parameters are shown in Figure 8. The results of theinversion using the line source reveal the heterogeneity of the faulting: (1) the
average rupture velocity is approximately 2.0 km/sec for the northern subfault and2.5 km/sec for the southern subfault, and (2) the relative moment release is highest
at the northern edge (element 7) on the northern subfault and at the southern edge
(element 1) on the southern subfault. Moreover, it is noted that the rupture

1632

EIICHI FUKUYAMA AND KOJIRO IRIKURA
{SEC)

GO
50
40
30
20
10
0

RUPTURESTARTING TIME

i:it
RELATIVE MOMENTRELEASE

FIG. 7. T h e results of the inversion for the line source. Dotted lines are the starting mode] and solidlines are the final model.

propagated smoothly in the southern subfault. These results are assured to bereliable judging from the remarkable decrease of the standard deviations (Table 2).
A r e a s o u r c e case. We performed forward modeling for a finite area source modelusing the same procedure as for the line source to determine the initial model forthe inversion. Here again, we fixed parameters m~ to 1 and controlled parameter t,
by varying both the northern subfault rupture velocity and southern subfault rupture
velocity. We assumed circular rupture propagation away from an initiation point.
From the results of forward modeling, we find the rupture velocities are again
approximately 2.0 km/sec along the northern subfault and 2.5 km/sec along thesouthern subfault, with the rupture initiating at element 2b (Figure 5b), giving
similar results as the inversion of the line source case. Then, we use the parameters
associated with these rupture velocities to construct the starting model of theinversion (Table 3). We assume the covariance matrices E and D defined by (A6)and (A7) are as follows: all off-diagonal elements of both E and D equal to zero. All
diagonal elements of E are equal to 0.3. Diagonal elements of D corresponding to
the relative moment releases of northern and southern subfaults and rupture
starting times of northern and southern subfaults are equal to 4.0, 2.0, 8.0, and 4.0,
respectively.
After approximately five iteration steps, we obtained the final solution, which
satisfies the convergence criterion of the inversion [ ( H A ) t r / ( K I ) t r = 0.05]. The
parameters of the starting model and final model are shown in Table 3, and thefinal results are shown in Figure 9. The waveforms which were calculated from thefinal model are shown in Figure 10. Since standard deviations of the area source
parameters are much greater than those for the line source parameters as shown in
Table 3, the results obtained here involve some uncertainties. This fact may bepartially caused by the increase of the parameter. Since the synthetic seismograms

RUPTURE PROCESS OF THE 1983 JAPAN SEA EARTHQUAKE

LINE

SOURCE

1633

MODEL

0.763

0.610

QIK N

OFU N

0.6|6

0.668
FIIK gOFU[v
O. 552

0.565

0.331

0.633

0.585

0.721

SEN N

....

~

'

"

• ~"
O, 686

O. 529
5£N g

HnK

q

15sec

U

E

J

0.822

suT .

0.922

fiUT C

FIG. 8. The waveforms calculated from the final parameters of the line source. Solid lines are observed
seismograms and dotted lines are synthesized. The upper left number of each trace is the value of thecross -correlation.

for the final model show a better agreement with the observed ones, we consider
these results to give the essential characteristics of the two-dimensional source
process. The inhomogeneous slip distribution and the nonuniform rupture propagation for the area source model are consistent with those of the line source model.
DISCUSSION

The approach we propose here is a first attempt in applying empirical Green's
functions to study fault heterogeneity. There are several aspects that need to berefined in our algorithm, i.e., the waveform synthesis and the inversion procedure.
Empirical Green's [unctions. We need to evaluate the validity and applicabilityof the empirical Green's function method. Criteria are needed to determine which
event is the most suitable to the empirical Green's function for each element on thefault surface. Here, we used the records from only two large aftershocks, assigningthem to the northern subfault and the southern subfault of the main shock surface,

1634

EIICHI FUKUYAMA AND KOJIRO IRIKURA

P"

C:)

¢¢

"~

¢D

"¢~

c'q

CD

5"q

~'-

t"'*

~'~

I ~

~~ . ~v-d I ¢0¢~r-d r-d

~

[ ~

r--4

I ~'~

~

o0

°5xl~ r-d I ~r-d r-~

~

c',l

c¢~

I~

~

i oo:~ ¸¢,D°~Z'~~ *-d I ~~

~-d

©

r-d

~

"~0

~'J

,<

O

O
~J
¢~

<DO

[-.-

~,~

¢q ~

~,,~

c.o ¢~

~

-~

~-- ~

~

,~

¢'q ~

C:~'~
¢D

I

-b-

RUPTURE PROCESS OF THE 1983 JAPAN SEA EARTHQUAKE

1635

NIO~W
"

"

4 6 sec(a) RUPTURE STARTING TIME

(m)
7.0
5.0

1.0~

3.0
0.0

(c) SLIP DISTRIBUTION
Fla. 9. The results of the inversion for the area source. (a) The rupture starting time on each element.
(b) Distribution of the relative moment release. (c) If we assume the total moment equals to 4.55 x 1027
dyne.cm {after Dziewonski et al., 1983) and rigidity equals to 3.5 x 1011 dyne/cm ~, we can translate (b)into the distribution of the slip.

respectively. We used the same records for all elements in each subfault as empiricalGreen's functions, correcting them only for the time delay that corresponds to therupture propagation over the fault surface and the seismic propagation from eachelement to the station, and for the amplitude due to the geometrical spreading
factor. This means that we assigned an average Green's function to each elementon the fault. The capability to resolve the irregularities on the fault surface isrestricted by the size of the available aftershocks. The aftershock of M = 7.1, which
was used as an empirical Green's function for the northern subfault, may be toolarge to resolve the fault motion in detail. To improve the resolution of the rupture

1636

EIICHI FUKUYAMA A N D KOJIRO IRIKURA

AREA

SOURCE

MODEL
0.735

0.79l
RIK N

"

0.719
I::IIK E. -,. . . .
0k.494 '

I

OFU N

0.594
i % ~,

"

i', n

r

/

~'

OFU E

0.558

.

0,679

O. 538

0,751

0.525
HRK N

SFN N

0,668

0.671
(t

sEN E

L ~ J

15sec

O. 800

5UT N ~ , . . ~ _ ~O. P025u,

Fro. 10. The waveforms calculated from the final parameters for the area source model.

process, we need to reconstruct the empirical Green's functions using some smaller
events that are distributed uniformly within the fault region.
Synthesis method. We can correct for amplitude variation caused by the radiation
pattern if we assume the wave type of the records and the focal mechanism. We
infer from the travel-time distance curve that the wave motions we treated in thisstudy consist of S waves as well as surface waves locally converted from S waves.
We evaluated the influence of radiation pattern on the amplitude variation by
calculating the ratio between the radiation pattern of each element and that of theaftershock, assuming that the observed seismograms contain mostly S waves. We
found that for almost all of the stations the ratios are within the range of 0.9 and1.1. Therefore, we can neglect the correction for this effect.
Since the waves of interest are mainly body waves (including surface waves locallyconverted from body waves) and reflected waves from the reflector inside the earth,p in equation (3), for correcting geometrical spreading factor, might be between 0.5
and 1.0. We calculate the terms (re/ri) and (re~r,) ~n for the worst case of our model,

RUPTURE PROCESS OF THE 1983 JAPAN SEA EARTHQUAKE

1637

and we find the difference is within the range of several per cents. This order is as
much as or less than that of the estimated errors of the radiation correction term.
Therefore, we choose the value o f p for the most remarkable phase, i.e., direct bodywaves.
Total moment. If we exactly know the seismic moments for the aftershocks usedas empirical Green's functions, we can calculate the total moment of the main shockusing equation (15).
Mo t°tal -~ m s~I N sJ~3Mdo61 + m n ( N n ) 3 M o 7"~(15)

where ms and m, are the means of m, which belong to the southern and northernsubfaults, respectively. Ns and Nn are the numbers determined by equation (2) for
the southern and the northern subfaults, respectively. However, in this inversion,
we do not estimate the absolute moment of the main shock but only the moment
ratio between main shock and its aftershock. According to Dziewonski et al. (1983),
if we adopt M071 = 1.9 × 102~ dyne. cm and M0 ~'1 = 9.6 × 1024 dyne. cm, we find thetotal moment of the main shock to be 3.0 x 1027 dyne. cm. This value is less thanthat listed in Dziewonski et al. (1983) but is comparable to that of Ishikawa et al.
(1984).
Model parameterization. We evaluate the inhomogeneity of the faulting process
by parameter m,, relative moment release on each subfault. We find that thisparameter is very sensitive to the waveform. The variation of parameter m, corresponds to the variation of the slip velocity. This means physically a variation of thestress drop. The rise time determined by equation (1) is fixed in this analysis
because the waveform used in the synthesis has predominant frequencies in therange lower than 0.2 Hz, where the value of the rise time has less influence.
Therefore, the variation of parameter mi may come from variations in the slip
displacement distribution.
Data form during inversion. The cross-correlation error function defined byequation (6) and used in this inversion procedure can be expressed after somemanipulationsify, = ~[Oi(t) = S,(t + t ~ ) ] 2 dt

whereC),(t) = o,(t)/[f o,(t), at]1''
S,(t) =

s,(t)/[fSi(t) 2 dt] 1/2

and Oi(t) and S,(t) are observed and synthesized seismograms of ith station,
respectively. This relation shows that minimizing the sum of squares of yi leads tominimizing the difference between normalized synthetic seismograms and normalized observed seismograms. Therefore, the information on the absolute amplitude
as well as the absolute time on each record is automatically removed, and we treat
only the waveform information. In our present analysis, it is proper to use the datain the form of cross-correlation error function because it is not so easy to calibrate

1638

EIICHI FUKUYAMA AND KOJIRO IRIKURA

the absolute time and absolute amplitude of each seismogram exactly. However,this may be an overreduction of the information contained in the data, especiallyfor the area source model. If we can use high quality data with the absolute timeand amplitude, we can extract more information using a different inversion procedure.
Results from the inversion. The bending corner of the fault is a very interestingregion. This region, located at the intersection of the Sado ridge and Okushiri ridge(Figure 4), may represent geologically inhomogeneous aspects. Nakamura (1983)
suggests the possibility of the existence of a young subduction zone between the
North American plate and the Eurasian plate, and this plate boundary crosses thisfault at the bending corner. From the results of this analysis (Figures 7 and 9), we
find that less moment was released, and rupture propagation was temporarily healedin this region during the main shock. This region may have acted as a geometricalbarrier or an inhomogeneous barrier (Aki, 1979) during the main shock.
CONCLUSIONS

The rupture process of the 1983 Japan Sea earthquake has been studied using anonlinear inversion based on the Bayesian approach by comparing the observed
seismograms with synthetic ones derived from empirical Green's functions. The
synthetic seismograms and the observed seismograms show a good agreement foralmost all of the stations after approximately five iterations. The moment releaseis maximum at the northern edge and the southern edge of the main shock fault,
and minimum around the intermediate region between the southern subfault andthe northern subfault of the fault surface. The average rupture velocity is approximately 2.5 km/sec for the southern subfault and 2.0 km/sec for the northernsubfault.
ACKNOWLEDGMENTS
The authors wish to thank Professors Soji Yoshikawa, Takeshi Mikumo, ¥oshimasa Kobayashi, andEgill Hauksson for constructive reviews of this manuscript. Professors Keiiti Aki and David D. Jacksonalso gave helpful advice and criticism. Data were provided by Japan Meteorological Agency, and Dr.
Shigeru Noda permitted us to use his digitized data of the SENDAI station.
This work was performed while Dr. Irikura was on leave from the Disaster Prevention ResearchInstitute, Kyoto University, as a visiting Professor at the Department of Geological Sciences, Universityof Southern California. It was partially supported by the National Science Foundation under Grant
CEE-8408227.

REFERENCES
Aki, K. (1967). Scaling law of seismic spectrum, J. Geophys. Res. 72, 1217-1231.
Aki, K. (1979). Characterization of barriers on an earthquake fault, J. Geophys. Res. 84, 6140-6148.
Dziewonski, A. M., J. E. Franzen, and J. H. Woodhouse (1983). Centroid-moment tensor solutions forApril-June, 1983, Phys. Earth Planet. Interiors 33,243-249.
Geller, R. J. {1976). Scaling relations for earthquake source parameters and magnitudes, Bull. Seism.
Soc. Am. 66, 1501-1523.
Hadley, D. M. and D. V. Helmberger (1980). Simulation of strong ground motions, Bull. Seism. Soc. Am.
70, 617-630.
Hartzell, S. H. (1978). Earthquake aftershocks as Green's functions, J. Geophys. Res. Letters 5, 1-4.
Hartzell, S. H. and T. H. Heaton (1983). Inversion of strong ground motion and teleseismic waveform
data for the fault rupture history of the 1979 Imperial Valley, California, earthquake, Bull. Seism.
Soc. Am. 73, 1553 1583.
Heaton, T. H. and D. V. Helmberger (1979). Generalized ray models of the San Fernando earthquake,
Bull. Seism. Soc. Am. 69, 1311-1341.
Irikura, K. (1983). Semi-empirical estimation of strong ground motions during large earthquakes, Bull.
Disas. Prey Res. Inst., Kyoto Univ. 33, 63-104.

RUPTURE PROCESS OF THE 1983 JAPAN SEA EARTHQUAKE

1639

Irikura, K. and I. Muramatu (1982). Synthesis of strong ground motions from large earthquakes using
observed seismograms of small events, Proc. 3rd Inter. Earthq. Microzonation Conf., Seattle 1,447458.
Ishikawa, Y., M. Takeo, N. Hamada, M. Katsumata, K. Satake, K. Abe, M. Kikuchi, K. Sudo, M.
Takahashi~ S. Kashiwabara, and N. Mikami (1984). Source process of the 1983 Japan Sea earthquake,
Chikyu 6, 11-17 (in Japanese).
Jackson, D. D. (1979). The use of a priori data to resolve nonuniqueness in linear inversion, Geophys. J.
R. Astr. Soc. 57, 137-157.
Jackson, D. D. and M. Matsu'ura (1985). A Bayesian approach to nonlinear inversion, J. Geophys. Res.
90, 581-591.
Japan Meteorological Agency (JMA) (1984). Report on the Nihonkai-Chubu earthquake, 1983, Technical
Report of the Japan Meteorological Agency, No. 106 (in Japanese).
Kanamori, H. (1979). A semi-empirical approach to prediction of long-period ground motions from greatearthquakes, Bull. Seisrn. Soc. Am. 69, 1654-1670.
Kanamori, H. and D. L. Anderson (1975). Theoretical basis of some empirical relations in seismology,
Bull. Seism. Soc. Am. 65, 1073-1095.
Kikuchi, M. and H. Kanamori (1982). Inversion of complex body wave, Bull. Seism. Soc. Am. 72,491506.
Maritime Safety Agency Hydrographic Department (1984). Results of sea bottom survey in 1983 Japansea earthquake hypocentral region, Report of The Coordinating Committee Earthquake Prediction,
Vol. 31, 104-107 (in Japanese).
Mellman, G. R. (1980). A method of body-wave waveform inversion for the determination of earth
structure, Geophys. J. R. Astr. Soc. 62, 481-504.
Mori, J. and K. Shimazaki (1983). Source process of the May 26, 1983 Japan sea earthquake, Abst. Meet.
Seism. Soc. Japan 2, A16.
Nakamura, K. (1983). Possible nascent trench along the Eastern Japan Sea as the convergent boundary
between Eurasian and North American plates, Bull. Earthquake Res. Inst., Tokyo Univ. 58, 711722 (in Japanese).
Ruff, L. and H. Kanamori (1983). The rupture process and asperity distribution of three great earthquakes
from long-period diffracted P-waves, Phys. Earth Planet. Interiors 31,202-230.
Sato, T. (1983). Source properties of the 1983 Japan sea earthquake as inferred from aftershock
distribution and strong motion accelerograms, Ph.D. Thesis, Tohoku University, Sendal, Japan.
Sato, T. (1985). Rupture characteristics of the 1983 Nihonkai-Chubu (Japan Sea) earthquake as inferredfrom strong motion accelerograms, J. Phys. Earth 33, 525-557.
Shimazaki, K. and J. Mori (1983). Focal mechanism of the May 26, 1983 Japan sea earthquake, Abst.
Meet. Seism. Soc. Japan 2, A15.
Spudich, P. and N. L. Frazer (1984). Use of ray theory to calculate high-frequency radiation fromearthquake sources having spatially variable rupture velocity and stress drop, Bull. Seism. Soc. Am.
74, 2061-2082.
Tarantola, A. and B. Valette {1982). Generalized nonlinear inverse problems solved using the leastsquares criterion, Rev. Geophys. Space Phys. 20, 219-232.
DISASTERPREVENTION RESEARCHINSTITUTE
KYOTO UNIVERSITY
UJI, KYOTO 611, JAPANManuscript received 19 February 1986
APPENDIX

We consider the inversion problem that data vector y can be expressed by themodel parameter vecotr x and random error vector e with the formy = f(x*)

+ e

(A1)

w h e r e f is t h e m o d e l f u n c t i o n a n d s u p e r s c r i p t * d e n o t e s t h e t r u e v a l u e . W e s o l v e(A1) u s i n g t h e B a y e s i a n a p p r o a c h , w h i c h w a s p r o p o s e d b y J a c k s o n a n d M a t s u ' u r a(1985). T h e B a y e s i a n a p p r o a c h h a s a n a d v a n t a g e t h a t a p r i o r i i n f o r m a t i o n , w h i c h1640

EIICHI FUKUYAMA AND KOJIRO IRIKURA

was found by other analysis, can be directly incorporated. A priori information isexpressed in the form of initial estimates of the parameters, Xo and the relationx o = x* + d(A2)

where d is an unknown error vector of a priori estimates. We incorporate a priori
information by Bayesian ruleP(xl y) = P(yi x).P(x)/P(y)(A3)

Here, P ( x l y ) is a posteriori joint probability density function for the parameters,
is equivalent to P(e), P ( x ) is a priori joint probability density function of
the parameters, and P ( y ) is a normalizing factor.
For the hypothetical solution x, we define the residuals of observational and apriori data as
P(yix)

e : y - f(x)

(A4)

d =

(A5)

x0 -

x.

Suppose these errors are Gaussian with zero means and covariance matrix E andD, respectively, i.e.,e - N(0, E)(A6)

d ~ N(0, D).

(A7)Then, a posteriori probability density function isP ( x l y ) = a exp(-T2/2)(A8)

T 2 = eTE-le + dTD-ld.

(Ag)

where

We can obtain the maximum likelihood solution by minimizing (A9). The solution
can be found by solving the implicit equation
A TE-le + D-ld = 0(A10)

A~ = (af~/ax~ L.

(A11)

where