Bulletin of the Seismological Society of America, Vol. 87, No. 4, pp. 904-917, August 1997

High-Frequency Radiation Process during Earthquake Faulting
Envelope Inversion of Acceleration Seismogramsfrom the 1993 Hokkaido-Nansei-Oki, Japan, Earthquakeb y Y a s u m a r o K a k e h i a n d K o j i r o Irikura

Abstract

We investigate the process of high-frequency (1 to 10 Hz) radiation on
the fault plane of the 1993 Hokkaido-Nansei-Oki, Japan, earthquake ( M w = 7.5)
from the envelope inversion of strong-motion acceleration seismograms. For the
analysis, empirical Green's functions are used because theoretical approach is notavailable for such high frequencies. The source is modeled with two fault planes
with different strike angles. The rupture process of this earthquake is very complex
in terms of high-frequency wave generation. The rupture, which started on the northern fault plane, had a delay of about 10 sec or propagated very slowly between thenorthern and southern fault planes. High-frequency wave radiation is large at thenorthern and southern edges of the source region. Deceleration of rupture is also
observed there. This is interpreted to be associated with stopping of rupture. Anotherhigh-frequency wave radiation area is found at the center of the northern fault plane,
where discontinuity in the depth distribution of aftershocks suggests an existence ofa barrier. The areas of high- and low-frequency wave radiation are not correlated.
This is considered to result from the complexity of rupture process. We cannot distinguish between westward and eastward dip of the southern fault plane because of
one-sided station distribution.

Introduction
plicated to explain with two or three fault planes. Nakanishiet aL (1993) obtained different mechanisms from an analysisof P-wave polarities and from a moment tensor inversion of
long-period (more than 70 sec) waveforms. The mechanism
from P-wave polarities showed a westward dipping fault
with a low angle (or eastward with a high angle), while theirlong-period wave analysis gave a mechanism dipping westward with a high angle (or eastward with a low angle), very
similar to the Harvard CMT solution. This difference is because the former reflects the initial stage of rupture, but thelatter describes the whole rupture process. Kuge et al. (1994)
analyzed both teleseismic body waves and long-period surface waves. From the analysis of the body waves, using
Kikuchi and Kanamori's (1991) method, they found that it
was a very complicated multiple shock, whose mechanism
changed during rupture propagation (Fig. 2). At the northernpart of the source region, the fault plane is dipping westward
with a low angle (or eastward with a high angle), while atthe southern part, it is dipping westward with a high angle(or eastward with a low angle). Their long-period surfacewave analyses also support this result. They found that different mechanisms were obtained depending on whether
they used for analysis both Rayleigh and Love waves or only

The Hokkaido-Nansei-Oki, Japan, earthquake ( M w =7.5) occurred at 13:17 GMT on 12 July 1993, very close to
Okushiri Island, southwest of Hokkaido, Japan. The earthquake, subsequent tsunamis, landslides, and fires resulted in
197 deaths and 42 missing people. As this earthquake occurred at the suspected North American-Eurasian Plateboundary, knowledge of its fault plane is important for understanding the tectonics of this region. However, it became
clear that the fault plane geometry was very complex anddifficult to identify because different fault models were suggested by different data, such as aftershock distribution, teleseismic body waves, surface waves, geodetic data, and tsunami data.
Hino et aL (1994) located aftershocks using OBS data
(Fig. 1). The aftershock epicentral distribution suggests thatthe source region is divided basically into two parts: thenorthern part striking SSW and the southern part striking
SSE. Figure 1 also shows the topography around the Okushift Island (Hydrographic Survey Team of Hokkaido-Nansei-Oki Earthquake, 1994) on which the aftershock epicenters are plotted. The topography has a change of strike at thesame place as the aftershock distribution. However, Hino etal. (1994) found that the hypocentral distribution is too com904

905

Y. Kakehi and K. Irikura

I

o

Lo

I

,<dip

t~#o
139 °

/

1400
~-3"

\

\,

42 °

42 °

o
o

I
139 °

140 °

S

!

i

139 °

1400

!

Depth

Figure 1. (Left) The aftershock hypocenters (21 July, 12:00 to 16 August, 23:59)
located from the OBS observation (Hino et al., 1994). Solid circles are OBS stations.
(Right) The topography around Okushiri Island (Hydrographic Survey Team of Hokkaido-Nansei-Oki Earthquake, 1994) with the aftershock epicenters by Hino et al.
(1994).

Rayleigh waves. They showed, from simple numerical simulations, that the difference was due to the mechanism
change during the source process. Umino and Hasegawa(1994) determined aftershock locations using the s P phaseand showed that their depth distribution dipped westward in
the southern region. Hashimoto et al. (1994) presented amodel composed of three fault planes (two for the southernpart), all of which are dipping westward, to explain geodetic
data on Okushiri Island. Their northernmost and central
planes have low dip angles (20° and 30 °) and the southernmost plane has a high dip angle (55°). From analysis of tsunami data, Imamura e t al. (1994) presented a model composed of two fault planes. Though they could not determine
the direction of dipping of the northern fault plane, theyreported that the southern fault plane must dip eastward inorder to explain the tidal record at Esashi, close to the southern source region.
The complexity of the rupture process is also easily seenin the strong ground motion acceleration seismograms observed within 200 km from the source region (Fig. 3). Sincenear-field data have better spatial resolution than teleseismic
data, analysis of these accelerograms is very appropriate for
estimating the complex rupture process of the HokkaidoNansei-Oki earthquake. Moreover, the determination ofhigh-frequency wave radiation areas on the fault plane willcontribute greatly to the understanding of the dynamics ofrupture, because high-frequency wave radiation is related tothe irregularity of rupture propagation.
Major developments in the 1980s in the source inversion techniques allowed derivation of detailed source processes from the analyses of seismic waveforms (e.g., Kikuchi and Fukao, 1985). In particular, the usage of near-fielddata made it possible to determine the rupture processes ofearthquakes with high accuracy and resolution (e.g., Hartzell
and Heaton, 1983; Fukuyama and Irikura, 1986; Mendoza
and Hartzell, 1989). However, the high-frequency coverageof source process studies by the analysis of seismic waveforms has been limited to about 1 Hz. To obtain further
understanding of the earthquake dynamic rupture process,analysis using higher-frequency waves is required. Also, for906

Y. Kakehi and K. Irikurathat are too oscillatory to be inverted. We applied thismethod to the 1993 Kushiro-Oki, Japan, earthquake (Mw =7.3) and obtained the interesting result that high-frequency(2 to 10 Hz) waves were radiated mainly from near the periphery of the fault plane, while the slip distribution (Takeo
et al., 1993), related to lower-frequency radiation, was concentrated in the center of the fault plane.
Here we apply this method to the 1993 Hokkaido-Nansei-Oki earthquake to estimate high-frequency wave radiation areas on the fault plane. We discuss the dynamics of the
rupture in terms of the results.

Total

50 km

Method~f

OKS

A detailed explanation of the envelope inversion
method is in Kakehi and Irikura (1996). Here we give onlya brief explanation. We use the envelopes of acceleration
seismograms as source inversion data. The envelopes arecalculated as the running root mean square of acceleration
seismograms with a given time window length. This smoothing procedure allows an inversion for frequency greater than
1 Hz, but it reduces the inversion resolution. The model parameters are the weights and rupture times of the meshes
into which the fault plane is divided.
The synthetic waveform of the target event is calculated
from superposition of the small-event seismogram used as
an empirical Green's function, following Irikura (1986). In
order to satisfy the scaling law, N 3 small-event seismograms
are superposed, where N is the common ratio of fault length,fault width, and slip between the target and small events.
From the scaling law, the following relation is obtained:

1

3'0

go

9'0 seeMoment-Rate

Figure 2. The spatial and temporal distribution ofsubevents obtained from the teleseismic body-wave
inversion (Kuge et al., 1994).
earthquake engineering, we need to know from what part of
the fault plane seismic waves higher than 1 Hz are radiated,because those frequencies have a strong effect on common
low-rise structures. Kakehi and Irikura (1996) developed a
new method to determine high-frequency wave radiation
areas on the fault plane from the envelope inversion of acceleration seismograms. In this method, we adopted the empirical approach for calculating the Green' s function because
it was impossible to obtain sufficient accuracy theoretically.
Envelopes are used as inversion data instead of waveforms

sur

o~
'7,

Here, 0 is the flat level of the displacement spectrum at low
frequencies and Mo is the seismic moment, c is the stressdrop ratio between large and small events:

Idili-'=.i~,.litla, i.,i ~ r p l v - r , , , , q ~..... ' . ..........
...... "~~(2)

c = AaL/Aas.

Subscripts L and S indicate large and small events. Further,z; . A KsN(1)

ULIUs = MoLIMos = c N 3.

==. SApNS~m............ dlli~........., , , ~" " - ~','~,.,.,.~mln~! F - , " ' :NS......... ]lillil&i,=~.~_tt,

ii,¢tl,i,i~-~...... ~ l ? l l ~ r r " : i ' " " ' t......... - ....

l

U ....................
Di.... ~"i.......
i

i

-

'
i

i

i

U

......

: "
t

i

i

D
i

=

i

J

i

I £ .....

ii................~i~.l,~.,~...........
~,.,.ml~,.inl~]lgll.,.

I

i

=

],

i

•

.

i

i

U . .,~,.a_~.~6~=.-.~,.=........ ~'r..-itl~

.

D
i

~

..... " ~ " " n ~=it'~'~''`''''ii.........
i

i

"
i

i

i

0 10 20 30 40 50 60 70 80 90100

0 10 20 30 40 50 60 70 80 90100

0 10 20 30 40 50 60 70 80 90100

sec

sec

see

Figure 3. The observed acceleration waveforms of the 1993 Hokkaido-Nansei-Oki
earthquake. They are bandpass filtered between 1 and 10 Hz.

907

High-Frequency Radiation Process during Earthquake Faulting--Envelope Inversion of Acceleration Seismograms
Table lP a r a m e t e r s o f the A f t e r s h o c k s U s e d as Empirically G r e e n ' s F u n c t i o n
Event
No.

Origin

1
2
3
4

Time

Lat.

Long.

(deg)

(deg)

Mj~ea~

M0(*1017N-m)

fc (Hz)L (kin)*Area

Station

43.10
42.72
43.17
41.96

139.23139.33139.01139.89

6.0
6.1
5.2
6.3

4.2
8.9
0.63
20

0.33
0.33
0.83
0.33

7.6
7.6
3.0
7.6

northern area
central area
northern area
southern area

HAK, SAP
SUT, HAK, SAP
SUT
SUT, HAK, SAP(GMT)

12 July 1994 14:45
12 July 1994 16:01
14 July 1994 17:38
7 Aug. 1994 19:42*L is a rough estimation of source size calculated by L = (nr~) 1/2, where r is obtained from Brune's (1970, 1971) formula: r = 0.373fl/f(a)northern area
= SUT,- H A K

O~

: ] E

=o,J"

~. S A P......

.II~.,I~......................

lv
.....

~

N

U ......
D ~

""

....

U

.w

D

-

0 10 20 30 40 50 60 70 80 90100
sec. . . . . . . . . ;a.~. . . . . . . .
:~.......

~......

:":S.....~ 'r*1]~w]~"lr.........

'. . . . . .

U.......... ,[l'~t~..a,...., .........

.....

1'o ~o 3'o 4'0 ~'o 6'0 7'o 8'o 9'OlOO

0 10 20 30 40 50 60 70 80 90100
sec

seo

central area(b)":. S U T

= HAK

L

=

"*

N

.

U

~

~..

D

-

i

i

~
,

=

t

i

,

i

=

i

t

0 10 20 30 40 50 60 70 80 90100
sec

(c)

o SAP

....da .........
N:::.,;:)dlJ~ . . . . . .

U....... ~D........

i

i

~i,

....

" r
i

-~o

i

. . . . . ~ILI~,u.,L . . . . . . . . . . ,~ . . . . . .

U . . . . . ~...~.~ . . . . . . . . . . . . . . . . . . . . . . . .

"~ . . . . . . . .

,

N

D

,

=

,

i

i

i

i

0 10 20 30 40 50 60 70 80 90100
sec.............
iJr~,, ....... ,...........

=

i

[

i

i

i

i

t

0 10 20 30 40 50 60 70 80 90100
sec

southern area~" S U T
co
"'rr

•

N

.

.

.

.......
.

.

s -..?~,.,.-:: I,,

, ~

,

,

"

'~w'rr',.r ......

°

"

' s

~

T

~:l'a

U .......
D~. S A P

~ HAK

(~

.......

,,~....

"/," ~ ~ ~ , ~ , , ~U"

DlO 2'o 3'0 io s'o 6'0 7o do 9'o~ooseeFigure 4.

i.... ~P'
r

i

J

"
i

....
~

t

J

'

i

i

0 102030405060708090100
sec

i

i

i

i

J

t

i

i

J

i

0 102030405060708090100
sec

T h e aftershock acceleration w a v e f o r m s u s e d as empirical G r e e n ' s f u n c tions. T h e y are b a n d p a s s filtered b e t w e e n 1 and 10 Hz. Events 3 (for SUT) a n d 1 (forHAK and SAP) are u s e d for the (a) northern area, event 2 is u s e d for the (b) centralarea, a n d event 4 is u s e d for the (c) southern area.

,

i

908

Y. Kakehi and K. Irikura

if their spectra follow the o9- 2 model, the following relation
is obtained:fi,L/As = cN.

(3)

Here A is the flat level of the acceleration spectrum at highfrequencies. From (1) and (3), we can obtain N and c.
Since the ratios in fault length, fault width, and slip
between the large and small events are not exactly the same,
we use NX × NW × N T instead of N 3. The fault plane of
the large event is divided into NX x NW subfaults, with NX
divisions along the strike direction and NW along the dip
direction. The size of the subfault corresponds to that of thesmall-event fault plane. In each subfault N T small-event seismograms are summed using the transfer function( N T - 1)n,Fo(t) = 6 ( t -tij) + lln'~

radiation intensity. Heterogeneity of the high-frequency
wave radiation is expressed by the distribution of w 0.
We correct for attenuation on synthetic accelerations.
Hoshiba (1993) obtained Q~-l, Q:-l, and Q;-I = Q~-I +Q~-I from his multiple-lapse time window analysis of earthquakes in Japan whose hypocentral distances and focal
depths are less than 120 and 40 km, respectively. His resultshows that Qt has a linear dependence on frequency for > 1
Hz. We u s e Qt = 100ffrom his result for the Aomori station,
which is closest to the source region of the Hokkaido-Nansei-Oki earthquake. The inversion is insensitive to the valueof Qt, as is verified by trying the case of Qt = 200f also.
Correction of the Q effect is performed as follows. First, theacceleration data are Fourier transformed by FFT. In the frequency domain, the Fourier amplitude is corrected considering the attenuation by Q effect. Then the corrected time
domain data are obtained by inverse Fourier transformation,

6[t

k=l

-t o - {(kto = (r o -1 ) r L I / { ( N T - 1)n'}],ro)/fl + So.

(4)

Data and Analysis(5)

We use the acceleration records from the JMA (Japan
Meteorological Agency) stations: Suttsu (SUT), Hakodate(HAK), and Sapporo (SAP). Four aftershocks are adopted as
empirical Green's function events. As shown in Table 1,
three empirical Green's functions, each of which covers thenorthern part, central part, and southern part of the source
area, are used for producing synthetic waveforms at eachstation. The acceleration records are bandpass filtered between 1 and 10 Hz for the analysis. Figure 4 shows theaftershock records used as empirical Green's function. Since
HAK did not record the mainshock, we use the acceleration
record from the station of BRI in Hakodate (Building Research Institute), which is 2.3 km from HAK of JMA, for themainshock. Figure 5 shows the spectral ratio between the

The synthetic waveform is obtained from superposition:
NX NW

A(t) = c. ~ , ~wo'(r/r o)'FO(t)*a(t).

(6)

i=1 j = lHere, zL is the rise time of the large event and fl is the Swave velocity, r 0, r, and r O are the hypocentral distances ofthe small event, the large event, and the subfault (i,j), respectively, n' is a parameter chosen to be much larger thanNT in order to remove the artificial period that appears in
relation to zL/(NT - 1). The asterisk denotes convolution.
w 0 and s O are the weight and rupture time of subfault (i,j),respectively. These are the model parameters to be obtained
from the envelope inversion. A non-negative constraint isgiven to w 0 using the algorithm of Lawson and Hanson(1974), because sign is meaningless due to the procedure to
take running root mean squares.
This inversion is nonlinear because the time variable s 0
is included in the model parameters and the calculation ofroot mean squares is nonlinear. Therefore, the inversion islinearized and solved iteratively. First, synthetic acceleration
waveforms are calculated by the empirical Green's function
method based on the model obtained in the previous stageof iteration. Next, envelopes of the synthetic waveforms arecalculated. Then w 0 and s 0 are determined to minimize thedifference between observed and synthetic envelopes in theleast-squares sense with a non-negative constraint for w 0.
The obtained w 0 and s 0 are regarded as the final solution onconvergence; otherwise, they are used as the initial modelfor the next iteration stage. As we are analyzing high-frequency accelerations, wij represents the weight of contribution of each mesh to the high-frequency acceleration envelope amplitude. Therefore, w 0 is called the acceleration

10

=
D

m

m

1.0

=
,,,.,,/i i illrll

0.1
0.1

\

I r lllwll

1.0

i i illllll
10

10 2

frequency(Hz)
Figure 5. The spectral ratio between the records ofevent 4 observed at the Hakodate station of JMA andat the one of BRI. Spectrum is calculated from vectoral summation of spectra for the three components.
Solid line marks a level of 0.7.

909High-Frequency Radiation Process during Earthquake Faulting--Envelope Inversion of Acceleration Seismograms

7

I,-100km° "l

~

~

I~

record of the number 4 event from the station of JMA andthe one from the station of BR1 at Hakodate. The spectrallevel between 1 and 10 Hz at the station of JMA is 0.7 times
as large as that of the station of BRI. This is considered toexpress difference of amplifications between them. We use
the mainshock record from the Hakodate station of BRI witha correction factor of 0.7.
Results from all kinds of data, such as aftershock distribution, mechanisms obtained from waveform inversion,
tsunamis, and geodetic data, agree in that the source is very
complex. Moreover, different kinds of data suggest different
fault plane models, particularly for the southern half of thesource area. However, we want to use as simple a model as
possible for source inversion. We use a simple model based
mainly on the aftershock distribution (Hino etaL, 1994) andthe result from multiple-shock analysis using body waves(Kuge et al., 1994) and abandon attempts for a model satisfying all the data. We adopt a model composed of two fault
planes because the source area is roughly divided into twoparts: the northern part striking SSW and the southern part
striking SSE, according to the epicentral distribution of aftershocks. For the northern part, all the data suggest a fault
plane dipping westward with a low angle. We determine astrike angle of 200 ° from the aftershock distribution, consistent with those of the first and second subevents of Kuge etal. (1994). Based on the dip angles (34 ° and 28 °) of the firstand second subevents of Kuge et al. (1994), we adopt a dip
angle of 30 ° .
For the southern part, we cannot find a conclusive evidence to determine whether it is dipping westward or eastward. Hashimoto et aL (1994) presented a model composedof three fault planes (two for the southern part), all of which
are dipping westward, to explain the geodetic data on Okushiri Island. On the other hand, the model presented by Imamura et al. (1994), from the analysis of tsunami data, hasa southern fault dipping eastward. Therefore, we try two

100km

"1

Figure 6. The station locations (triangles)
and the two different fault plane models usedin this study. The model on the left has thesouthern fault dipping westward with 30° dip
angle; the model on the right has the southernfault dipping eastward with 55° dip angle.
Solid star shows the rupture starting point.
Open stars are the epicenters of the empiricalGreen's function events.

different models for the southern fault; one dipping westward and the other eastward. Dip angles of 55 ° (strike angleof 160 °) and 30 ° (strike angle of 340 °) are adopted, because
Kuge et aL's (1994) result suggests a large westward dip
angle and small eastward one. Though their westward dip

event 3(for SUT)
event i(for HAK, SAP)

1event 2for SLIT, HAK, SAP)

event 4for S UT, HAK, SPur~)

Figure 7. The three areas of the fault plane(marked with different tones), in each of which different empirical Green's function event is assigned.

910

Y. Kakehi and K. Irikura

angles for the third and fourth subevents are 73 ° and 72 °, we
use a smaller value of 55 ° consistent with the depth extent
of the aftershock distribution of Hino et al. (1994). The fault
plane sizes are determined based on the epicentral distribution and depth extent of aftershocks as follows: 80 × 42 km
for the northern plane, 72 × 48 km for the southern planedipping westward, and 72 × 42 km for the one dipping
eastward. Our fault model is shown in Figure 6. The depthof the rupture starting point, denoted by a solid star, is 23
km. Our two-fault-plane models are intentionally simple andare not necessarily consistent with Hino et al.'s (1994) interpretation of aftershock distribution or Kuge et al. 's (1994)
result. As mentioned later, we cannot determine whichmodel is better. This may be because the stations are locatedon the same side of the source region, and sufficient resolution is not obtained for the direction of fault width.
Table 1 lists the aftershocks used for the analysis, theirseismic moments, and rupture durations. These seismic moments and durations are obtained from the comparison between the observed and synthetic displacements. The tablealso lists a rough estimation of the source sizes, as calculated
from the comer frequencies (reciprocals of durations) using
Brune's (1970, 1971) formula. As explained in the previous

Obs.
Syn.
SUT

i

~

i

,

~

i

i

i

HAK

J

J

i

~

i

~

section, when the target event is a single shock, the parameters N and c used for making synthetic waveforms fromsmall-event records are determined from the moment ratio
and high-frequency spectral ratio. The Hokkaido-NanseiOki earthquake is a very complex multiple shock, and although we need to separate it into subevents and calculate
N and c for each of them, that is not possible because theseismic waves from different subevents are not separated but
superposed on the near-field seismograms. Therefore, wetake a different approach. First, NX × NW, numbers of division into subfaults, are determined 10 × 6 for the northernfault plane and 8 × 6 for the southern one based on therough estimation of the aftershock source sizes in Table 1.
In Figure 7, the three areas covered by the different empiricalGreen's function events are shown. The northernmost area
is covered by event 3 for SUT and by event 1 for HAK and
SAP. As event 3 is about half the size of event l, in order to
use the same subfault size, we prepare synthetic waveforms
by summing up 2 × 2 × 2 records of event 3 based on thesame procedure and use them as the Green's function. Next,
NT and c are determined as follows. As mentioned above, it
is impossible to take moment and high-frequency spectral
ratios for each subevent; it is possible only for the first sub-

~

i

i

~

i

,

i

SAP

,

i

i

i

i

i

p

t

,

~

i

i

i

~

i

i

0 10 2 0 3 0 4 0 5 0 6 0 7 0 8 0 9 0 1 0 0

0 102030405060708090100

0 10 2 0 3 0 4 0 5 0 6 0 7 0 8 0 9 0 1 0 0

sec

sec

sec

Obs.
Syn. _ _
SUT

i].w

l.w

:l.w

HAK

SAP

a~.

~I

UD

. .20.
0. 10

NS

• UD

. 40
. .so. 60. 70. .80. 90100. . . . .0 . 10. 20. 30. .40. so 60 70 00 90100
sec

sec

0

1'0'20 3'0'0'0'
4 s

sec

60

Figure 8. (Upper) Comparison of the observed envelopes and the synthetic ones forthe initial model for the westward-dipping southern fault. (Lower) Comparison of theobserved envelopes and the synthetic ones for the final solution. The rms window length
used for the calculation of envelopes is 5 sec. Two peaks corresponding to the waves
from the northern and southern fault planes are indicated with arrows.

' 0'0 90100

High-Frequency Radiation Process during Earthquake Faulting--Envelope Inversion of Acceleration Seismograms~= SUT,,o HAK

Syn. ~]~
I

E

ql.I

W

o

~

~

" Syn.
,' '"r"~i'rlTlq~B''.ll

, . •

W

0 ~0 :~0 30 40 ~0 60 40 80 ~0~00''*"~'~W"'~"' '.:!,, .,,~,,

w

Syn.

Syn.

s

•

........

0 40 :/0 3'0 ,~0 ~'0 6'0 ~0 8'0 901~0
u. . . . . . . . . . . . . . . . . . . . .~. . . . . . . . . . . . .

,D .-~.~.r,,-,-...................

Syn.

s

¢~ 0 40 20 30 ,~0 5'0 6'0 ~0 8'0 90~00

sec

sec

-~

.~. . . .

O)

0~,
-~

O)

Syn.
1'0 2'0 30 4'0 ~'0 6'0 7'0 8'0 9'0100

sec

0 1'0 20 3'0 ,~0 s'o 6'o ~0 80 9'0100

secu ~.~,,~al~,=,,.~..,~.,,.......

D "~"P~rRII~I~'z"I"qF''. . . . . Iv"...

Syn.

D ....'r~"~r . . . . . . . . .~" . C. . . . . . . . . . . .

40 20 3'0 4'0 s'o 6'0 ~0 8'0 9'0100~".'~l,r~.~",,,','rl."~0 60 40 ~0 90~00
O~ 0 ~0 20 :~0 40 sec: Obs. I I,.hll.
i t ~ ~

sec"~ N .
-N . ..~..,.=u~,.a=m=,m,,~..........
i]]1~ sec ::':":"" 'it Obs.
l,.,|,~J,~=IlL~~,~o

SAP~ Syn. I~'~o 0 40 :~0 :~0 40 50 60 70 eo 90~0

Syn. ~,~ ,LI$,,J,

911

D"'TI~II~I"P r ~ ' ' r ~ ' q ' ~1'0 2'0 3'0 4'0 5'0 6'0 7'0 8'0 9'0100

seo

sec

Figure 9. Comparison of the observed acceleration waveforms and the synthetic
ones for the final solution for the westward-dipping southern fault.

event because seismic waves from it are little contaminatedby those from subsequent subevents even in the near-fieldseismograms. W e calculate the moment ratio between thefirst subevent of Kuge et al. (1994) and our event 2 for thearea around the rupture starting point and obtained a valueof 110. The spectral ratio is calculated for the first 20 sec of
S waves of the mainshock and event 2 at each station. The
spectrum o f each event is calculated from vectoral summation of spectra for the three components. The average ratioat 1 to 10 Hz is 9.1. From (1) and (3), N - 3.5 and c - 2.6
are obtained. In the actual simulation, after several trials, we
use NT = 4 and c = 2.3 for the part covered by event 2.
After fixing c for event 2, we determine c for other events
as follows. For event 1, we select a subfault at the boundary
between the northern and central areas. Then we calculate
the contributions of that subfault to the records of SAP and
HAK for two cases: one with event 1 as the empirical Green' s
function and the other with event 2. The value of c for event
1 is determined so that the amplitudes at SAP and SUT are
the same for both cases. In this way, we obtain c = 0.71,
0.23, and 1.2 for events 1, 3, and 4, respectively. W e use NT= 4 for events 1 and 4, and for the Green's functions madefrom event 3. Other parameters used for the analysis are TL= 3.2 sec, fl = 3.7 km/sec, and n' = 60.
Fukuyama (1994) estimated the source-time function of
event 4 by deconvolution using smaller-event records. Theirresult is that event 4 is a multiple shock composed of two
small and one large subevents. Iwata et al. (1994), however,from the waveform inversion of near-field displacements,
reported that it is composed of two subevents, one is muchsmaller than the other. They also simulated the near-field
acceleration waveforms by forward modeling using the empirical Green's function method and found that the accelerations were well explained by the main event only. Therefore, it is acceptable to take event 4 as a single shock, as far
as high-frequency accelerations are concerned.
In order to reduce the number of model parameters, we
treat four subfaults (subfanlt size corresponds to that of theempirical Green's function event) as one mesh for envelope
inversion as shown in Figure 7. The initial model of acceleration radiation intensities is the uniform distribution overthe fault plane. The initial values of rupture times are forradial rupture propagation with a constant rupture velocityof 2.8 km/sec. The rms window length used for the calcu-

912

Y. Kakehi and K. Irikura

lation of envelopes is 5 sec. As the inversion is nonlinear as
mentioned above, it is solved iteratively.

of the observed and synthetic envelopes for the final solution
for the case of the southern fault dipping eastward. Its residual is 21%, where we take that for the initial model of theother fault model as 100%. Since the residuals of the twomodels are very similar, we cannot conclude which model
is to be preferred. This may be because the stations are located on the same side of the source region, and the dataset
does not have enough resolution for the direction of the fault
width. Figure 11 shows the final solutions for the two models. The acceleration radiation intensities are shown withtones, where darker tones indicate larger intensities. The distribution of acceleration radiation intensity is heterogeneous.
This means that, as was expected from the observed seismograms, the rupture process of this earthquake was complex in terms of high-frequency wave radiation. Large intensities are seen at the northern and southern edge of thesource region. We can interpret the high-frequency wave
radiation in these areas as associated with stopping of the
rupture caused by barriers. Moreover, the rupture times inthese areas are 1 to 4 sec larger than the initial values; thatis, the rupture propagation decelerated. This also suggeststhat the rupture stopped due to barriers.
Another area at the center of the northern fault or a little
to the north of the center also has a large high-frequency
radiation. Since it is not located at the edge of the source
region, we cannot relate it to a stopping of the rupture. But
this area is coincident with a clear discontinuity in depth
distribution of aftershocks (Hino et al., 1994) (Fig. 12). The
existence of this discontinuity, which suggests a strong heterogeneity in the medium, is reliable because they determined the aftershock hypocenters from the records of OBS
covering the source region. We propose that this heterogeneity is the cause of the high-frequency radiation. In thisarea in the northern half part of the northern fault plane, there
is not only large radiation of high-frequency waves but alsohigh aftershock activity. The may be associated, but the relation is not simple because there is an area with small highfrequency radiation and high aftershock activity on thesouthern fault plane.

Results and Discussion
As discussed above, we try two kinds of fault planemodels. First, we show the envelopes for the case of thesouthern fault dipping westward. In the upper panel of Figure 8, the observed envelopes and the synthetic envelopes
for the initial model are compared. Here, it should be notedthat the observed envelopes at SUT and HAK have two peaks
indicated with arrows, with the second one much larger than
the first at HAK. If we assume continuous rupture propagation between the northern and southern fault planes, we cannot explain the arrival time of the second peak. Therefore,
we put a delay of rupture between the northern and southernfault planes in the initial model. We tried inversions withdelays of 6, 8, 10, 12, and 14 sec and chose the delay of 10
sec, which brought the minimum residual. In the lower panelof Figure 8, the synthetic envelopes for the final solution are
shown with the observed. The shapes of the envelopes at thefirst peak at SUT are improved. The amplitudes of the second
peak at HAK increased to fit the observed. At SAP, the duration of the main bump, which is made of superposition ofwaves from the northern and southern fault planes, increased
to fit the observed. The residual is reduced to 20% of theone for the initial model. Here, residual is defined as
3Res. =

f {Oi~t(t) - Sist(t)}2/M2~tdt,~(7)

i st = 1

where Oist(t) and Sift(t) are the observed and synthetic envelopes at the ist-th station, respectively. Mist is the maximum amplitude of the Oist(t). In Figure 9, the synthetic acceleration waveforms for the final solution are comparedwith the observed. The features of the observed waveforms
are reproduced very well. Figure 10 shows the comparison

Obs.
Syn. _ _
SUT

HAK

=] uo

uo

O

~

i

O

i

i

i

i

i

i

i

i

i

i

0 10 20 30 40 50 60 70 80 90100
sec

SAP

J

i¸

i

K

i

i

L

i

i

i

I

i

0 10 20 30 40 50 60 70 80 90100
sec

L

t

i

i

,

i

i

~

i

i

,

0 10 20 30 40 50 60 70 80 90100
sec

Figure 10. Comparison of the observed envelopes and the synthetic ones for the
final solution for the eastward-dipping southern fault.

High-Frequency Radiation Process during Earthquake Faulting--Envelope Inversion of Acceleration Seismograms913

f

a c c e l e r a t i o n radiation intensity0~) 0,2 {).4 ~]6 o.g t,0rupture t i m e p e r t u r b a t i o n f r o m intial v a | u e
-1

0

+I

+2

+3 +4 (s)Figure 11. The final solutions for the westward-dipping southern fault (left) and forthe eastward-dipping one (right). For each solution, the left and right panels show theacceleration radiation intensities and rupture times (in seconds), respectively. Rupture
time values are not shown at the meshes with small radiation acceleration intensities
because the contributions of those meshes are small and the rupture times have low
reliability.

At the area where the fault changes the strike direction
from SSW to SSE, radiation of high-frequency waves issmall, though the change of strike direction is expected tobehave as a barrier for the rupture. In this area, as mentioned
above, we had to put a delay of rupture of 10 sec in theinitial model. There is a possibility that the rupture propagated very slowly in this area rather than paused for 10 sec,
resulting in small radiation of high-frequency waves.
Kakehi and Irikura (1996) found a complementary distribution of high- and low-frequency wave radiation on thefault plane of the 1993 Kushiro-Oki, Japan, earthquake (Mw= 7.3)..The high-frequency radiation was large near the
periphery, while the low-frequency radiation was large at
the center. This is explained by a simple crack model. Figure
12 shows the comparison of our result with the moment

release distribution by Japan Meterological Agency (1995).
They also used a model composed of two fault planes. The
strike and dip angles of their northern and southern fault
planes are (195 °, 35 °) and (165 °, 35°), respectively. The
lengths and widths of the two planes are both (75, 39 kin).
Their result is from the waveform inversion of strong grounddisplacements with predominant period of about 20 sec, sothat it can be interpreted as the distribution of low-frequency
wave radiation. There is no clear relation between the distributions of the high- and low-frequency radiation, different
from the case of the 1993 Kushiro-Oki earthquake. We interpret this as due to the strong complexity of the rupture.
We perform a resolution and stability test of the inversion as follows. First, we calculate the synthetic envelopes
using the final solution obtained above (we call this final

914

Y. Kakehi and K. Irikura

I

I

I

°

6

e

I

1 *

I

+1

i

ql

"%,t .,~+~ ,¢j, o
e

80

%,o

o

-

o

.,,

"4'..," •

I-,,

4~

.

h o

<3

o

I 1~39°" ' ~*II

5xl0~gNm

I&, go -~
o

o

q:~

O

Depth
Figure 12. Comparisonof the acceleration radiation intensity distribution for the
westward-dipping southern fault (left) and the aftershock distribution by Hino et al.
(1994) (middle). The clear discontinuity in depth distribution is indicated with an arrow.
The distribution of the seismic moment release (Japan Meterological Agency, 1995) isalso shown for comparison (right).
solution "true solution" hereafter). Then we add 10% random noise to these envelopes. By regarding these noiseadded envelopes as the observed data for the test, we perform an inversion to recover the true solution. We use thefault model with the southern fault dipping westward. The
initial model for the test inversion is the same as was used
for the actual earthquake records. We define the percentageof the noise as the percentage of the area between the noiseless and noise-added envelopes (see the upper panel of Fig.
13) to the area of the noiseless envelopes. The middle panelof Figure 13 shows the noiseless (solid lines) and noiseadded envelopes (thin lines). The synthetic envelopes for thesolution obtained from the inversion of the noise-added envelopes are shown with thin lines in the lower panel of Figure 13, where the noise-added envelopes are shown withsolid lines.

We perform six inversions using different random
noises. Figure 14 shows the mean values and the standard
deviations of the model parameters obtained in the six inversions. The values in the parentheses are the ones used in
the noiseless envelopes, that is, true solution. The resolution
and stability is judged by whether these values are recovered
well and stably. The meshes whose model parameters satisfy
both of the following criteria are drawn with solid lines: (1)the true solution in the parenthesis lies within "mean value+ / - standard deviation," and (2) the standard deviationsof the acceleration radiation intensity and rupture time perturbation is less than 0.1 and 0.3 sec, respectively.
The model parameters of the northern half part of thesouthern fault are not recovered well, but they are inverted
stably, and the acceleration radiation intensities are estimated to be small, which is consistent with the tendency of

High-Frequency Radiation Process during Earthquake Faulting Envelopelnversion of Acceleration Seismograms

915

amplitudenoiseless

~

"

-

I. t

noiselessnoise-added
SUT

~TNS

o-I

HAK

~

~

¢~,

~
,

SAP

I~

o
i

i

i

,

i

=

l

i

,

,

0 10 20 30 40 50 60 70 80 90100

0 10 20 30 40 50 60 70 80 90100

Obs. _ _Syn. _ _
SUT

HAK

O

i

t

i

i

i

i

,

i

,

0 10 20 30 40 50 60 70 80 90100
sec

soo

see

t

NS

'1

0 10 20 30 40 50 60 70 80 90100
sec

SAP

4

i

i

,

i

,

i

~

i

i

i

,

0 10 20 30 40 50 60 70 80 90100
sec

,

i

i

L

i

i

i

i

,

i

i

0 10 20 30 40 50 60 70 80 90100
sec

Figure 13. The percentage of the noise added for the resolution and stability test of
the inversion is defined as that of the area between the noiseless and noise-addedenvelopes to the area of the noiseless envelopes (upper). The noiseless (solid lines) andnoise-added envelopes (thin lines). The percentage of the noise is 10% (middle). The
noise-added envelopes (solid lines) and the synthetic envelopes for the solution obtained
from the inversion of the noise-added envelopes (thin lines) (lower).

the true solution. Therefore, it is reliable that the high-frequency wave radiation is not large in this area.
Madariaga (1977) showed theoretically that high-frequency waves are radiated where the rupture velocitychanges abruptly. Stopping of rupture is an extreme case. In
this study, however, rupture velocity change is not seen in
the high-frequency radiation area coincident with the aftershock depth discontinuity. A possible reason is that the meshsize for the inversion is not small enough to detect such
rupture velocity change. The relatively large mesh size results from the small number of available stations.
We cannot yet discuss either general characteristics orregional differences about high-frequency wave radiation
because only two earthquakes have been analyzed using thismethod. Such a discussion requires an increase in the number of earthquakes analyzed. Interesting questions include,for example, whether there is a difference between interplateand intraplate earthquakes and whether there is a dependence
on magnitude.
By applying our method to high frequencies and the
waveform inversion method to low frequencies, we canmake a superbroadband (0.01 to 10 Hz, for example) analysis of an earthquake. This will enable us to reproduce thecomplete rupture history, which will contribute to the understanding of earthquake rupture dynamics.
Acknowledgments
Fujio Kusano of Japan Meteorological Agency and Yoshikazu Kitagawa
and Toshihide Kazima of Building Research Institute kindly provided us
with the acceleration records. Discussions with Tomotaka Iwata were helpful throughout this study. The manuscript was greatly improved by thevaluable comments of the two reviewers. Dr. Alan Linde kindly checkedand corrected the English. We are grateful to all of them. This research was
supported by the JSPS Fellowships for the Japanese Junior Scientists.

916

Y. Kakehi and K. IrikuraNE
0.83_+0.10(0.69)
0.80+0.06(0.82)

acceleration radiation intensity
0.39±0.02(0.39)
0.05_+0.06(0.00)

0.024-0.02(o.oo)
0.58±0.02(0.61)

0.31+0.04

0.92-+0,08(1.00)

0.03±0.02(0.00)
0.50+0,02(0.50)

0.69_+0.13

0.60_+0.20

0.20__+0.07

0.05-+0.04(1.00)(0.69)(0.10)(0.00)

0.50_+0,18(0.45)

SE(0.24)(0.32)(0.39)

O.50+O.O5(O.37)

0.25--0.18

0.21±0.08

0.17+__0.10

0.56±0.06(0.37)(0.10)(0.25)(0.52)

0.34_+0.05

0.33_+0.06

0.17_+0.06

0.63±0.12

1.00_+0.01(0.30)(0.48)(0.00)(0.69)(1.00)

0.35+0.12

0.24_+0.07

NW

SW

NE
0.0±0.3(o.o)

1.7±0.8(1.0)

1.1±0.3(t .0)

0.0-0,20.0)

0.4±0.5

0.9±0.9

(1.0)

(0.0)

rupture time perturbation (s)

SE

(-)

0.0+0,2(o.o)

0.0±0.3(o.o)(-)

0.0±0.0
0.0)

0.0±0.3(-1.0)--

0.0---0.1

(-)

(-)

(13.0)(-)

0.0±0.1(o.o)

0.0±0.0(0.0)

0.3±0.8

0.2±0.4

~.~

(-)

(2.0)
0.2±0.7(1.~

1.0±02(1.~

0.7±0.5-

0.0±0.3

3.0±0.3

(1.~

(-)

(0.~

(3.~

NW

SW
F i g u r e 14.
The mean values and standard deviations of the model parameters obtained from the inversions of the six different cases of the envelopes with 10% random
noise.

References
Brune, J. N. (1970). Tectonic stress and the spectra of seismic shear waves
from earthquakes, J. Geophys. Res. 75~ 4997-5009.
Brune, J. N. (1971). Correction, J. Geophys. Res. 76, 5002.
Fukuyama, E. (1994). Source processes of the 1993 Hokkaido-Nansei-Oki
earthquake and its largest aftershock, Kaiyo Monthly, Special Edition
for the Hokkaido-Nansei-Oki Earthquake and Tsunamis, 29-34 (inJapanese).
Fukuyama, E. and K. Irikura (1986). Rupture process of the 1983 Japan
Sea (Akita-Oki) earthquake using a waveform inversion method, Bull.
Seism. Soc. Am. 76, 1623-1640.
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.
Hashimoto, M., S. Ozawa, A. Yoshimura, T. Sagiya, T. Tada, and H. Tsuji(1994). Crustal movements associated with the 1993 Southwestern
off Hokkaido earthquake and its fault model, Kaiyo Monthly, Special
Edition for the Hokkaido-Nansei-Oki Earthquake and Tsunamis,55-61 (in Japanese).
Hino, R., T. Kanazawa, K. Suyehiro, T. Sato, and H. Shimamura (1994).
Aftershock distribution of the 1993 southwest off Hokkaido earthquake by ocean bottom seismographic array observation, Kaiyo
Monthly, Special Edition for the Hokkaido-Nansei-Oki Earthquake
and Tsunamis, 35-42 (in Japanese).
Hoshiba, M. (1993). Separation of scattering attenuation and intrinsic absorption in Japan using the multiple lapse time window analysis offull seismogram envelope, J. Geophys. Res. 98, 15809-15824.
Hydrographic Survey Team of Hokkaido-Nansei-Oki Earthquake (1994).
Preliminary results of the survey on the source region of the Hokkaido-Nansei-Oki earthquake, Report of Hydrographic Res. 30,
395-412 (in Japanese).
Imamura, F., T. Takahashi, and T. Takahashi (1994). Is the fault planedipping toward the east or west?--interpretation from tsunami data,
Kaiyo Monthly, Special Edition for the Hokkaido-Nansei-Oki Earthquake and Tsunamis, 179-184 (in Japanese).

Irikura, K. (1986). Prediction of strong acceleration motion using empiricalGreen's function, Proc. 7th Japan Earthquake Eng. 8, 37-42.
Iwata, T., K. Kamae, and K. Irikura (1994). Source process of the largest
aftershock (1993/8/8 MjMA 6.3) of the 1993 Hokkaido-Nansei-Oki
earthquake estimated from the near field strong ground motions,Kaiyo Monthly, Special Edition for the Hokkaido-Nansei-Oki Earthquake and Tsunamis, 80-87 (in Japanese).
Japan Meteorological Agency (1995). Report on the Hokkaido-Nansei-Oki
earthquake, 1993. Technical Report of the Japan Meteorological
Agency 117 (in Japanese).
Kakehi, Y. and K. Irikura (1996). Estimation of high-frequency wave radiation areas on the fault plane by the envelope inversion of acceleration seismograms, Geophys. J. lnt. 125, 892-900.
Kuge, K., M. Kikuchi, and J. Zhang (1994). Complex source process of theJuly 12, 1993, Hokkaido Nansei-oki eatlhquake estimated from teleseismic body waves and surface waves, Kaiyo Monthly, Special Edition for the Hokkaido-Nansei-Oki Earthquake and Tsunamis, 21-28(in Japanese).
Kikuchi, M. and Y. Fukao (1985). Iterative deconvolution of complex bodywaves from great earthquakes--the Tokachi-Oki earthquake of 1968,
Phys. Earth Planet. Interiors 37, 235-248.
Kikuchi, M. and H. Kanamori (1991). Inversion of complex bodywaves--III, Bull. Seism. Soc. Am. 81, 2335-2350.
Lawson, C. L. and R. J. Hanson (1974). Solving Least Square Problems,
Prentice-Hall, Englewood Cliffs, New Jersey.
Madariaga, R. (1977). High frequency radiation from crack (stress drop)
models of earthquake faulting, Geophys. J. R. Astr. Soc. 51, 625-651.
Mendoza, C. and S. H. Hartzell (1989). Slip distribution of the 19 September 1985 Michoacan, Mexico, earthquake: near-source and teleseismic constraints, Bull. Seism. Soc. Am. 79, 655--669.
Nakanishi, I, S. Kodaira, R. Kobayashi, M. Kikuchi, and M. Kasahara(1993). The 1993 Japan Sea Earthquake: quake and tsunamis devastate small town, EOS 74, 377-379.
Takeo, M., S. Ide, and Y. Yoshida (1993). The 1993 Kushiro-Oki, Japan,
earthquake: a high stress-drop event in a subducting slab, Geophys.
Res. Lett. 20, 2607-2610.
Umino, N. and A. Hasegawa (1994). Aftershock distribution of 1993 south-

High-Frequency Radiation Process during Earthquake Faulting Envelope lnversion of Acceleration Seismograms
west off Hokkaido earthquake estimated from near field sP phase,

Kaiyo Monthly, Special Editionfor the Hokkaido-Nansei-OkiEarthquake and Tsunamis, 49-54 (in Japanese).
Department of Earth and Planetary Sciences
Faculty of Science
Kobe University
Rokkodaicho 1-1, Nada-ku
Kobe 657, Japan
Tel: 81-78-803-0564; fax: 81-78-803-0490
E-mail: kakehi@ sc 122.earth.s.kobe-u.ac.jp(Y.K.)

Disaster Prevention Research InstituteKyoto University
Gokasho, Uji, Kyoto 611, Japan
Tel: 81-774-38-4060; fax: 81-774-33-5866
E-mail: irikura@egmdpri01.dpri.kyoto-u.ac.jp(K.I.)
Manuscript received 11 November 1994.

917