Bulletin of the Seismological Society of America, Vol. 98, No. 1, pp. 180–197, February 2008, doi: 10.1785/0120060183

Strong-Motion Generation Areas of a Great Subduction-Zone
Earthquake: Waveform Inversion with Empirical Green’s
Functions for the 2003 Tokachi-oki Earthquake
by A. Nozu and K. Irikura

Abstract

A waveform inversion with empirical Green’s functions was conducted
to estimate strong-motion generation areas (SMGAs) of the 2003 Tokachi-oki, Japan,
earthquake (Mw 8.3). Although the rupture process of this great subduction-zone
earthquake has been revealed with waveform inversions based on low-frequency(lower than 0.2 Hz) ground motions, it is much more important from an engineering
point of view to investigate how strong ground motions with higher frequencies weregenerated from the earthquake. Waveform data with higher frequencies have been
excluded from the conventional inversions because of the difficulty in computing
realistic theoretical Green’s functions at higher frequencies. In this study, with theaid of empirical Green’s functions, we extended the applicability of the waveform
inversion to higher-frequency ground motions up to 1 Hz. We selected three aftershocks for use in the analysis, referring to the similarity of the group delay time between the mainshock and the aftershock records. The resultant slip model has three
SMGAs, each of which is close to the hypocenter of one of the aftershocks used. Generally speaking, locations of the SMGAs identified in this study agree well with asperities identified from the inverted results using low-frequency (lower than 0.2 Hz)
ground motions plus geodetic data and tsunamis. It implies that strong ground motions
up to 1 Hz were generated from almost the same asperities producing lower-frequency
ground motions. In terms of the complexity of slip models, although our analysis isfocused on high frequencies, our slip model is at least as simple as conventional lowfrequency slip models. Such results would be useful in constructing source models of
future great subduction-zone earthquakes for strong-motion prediction.

Introduction
et al., 2000). Strong-motion data from the earthquake are far
superior to those from any past M 8-class earthquakes interms of both quality and quantity. The 1968 Tokachi-oki
earthquake (MJ 7:9), the source region of which is different
from that of the 2003 event, for example, was observed at
approximately 15 strong-motion stations (e.g., Tsuchida
et al., 1969). The 2003 Tokachi-oki earthquake, in contrast,
was observed at more than 600 K-NET and KiK-net stations(Honda et al., 2004). This amount of strong-motion data
gave us a rare opportunity to investigate how strong groundmotions were generated from a great subduction-zone
earthquake.
Slip models for the earthquake have been proposed byseveral authors. Yamanaka and Kikuchi (2003) investigated
the source process of the earthquake using teleseismic bodywaves. Yagi (2004) conducted a joint inversion of teleseismic
body waves and strong ground motions. Honda et al. (2004)
conducted a waveform inversion of strong-motion data.

A great subduction-zone earthquake occurred on 26September 2003 off the Tokachi Region, Hokkaido, Japan. According to the Japan Meteorological Agency (JMA) unified
source catalog, the hypocentral parameters are origin time,
04:50:07.64 JST; epicenter, 41.780° N, 144.078° E; depth,
42 km; and magnitude MJ 8:0 (Mw 8:3). As reported by
Hatayama et al. (2004) and Koketsu et al. (2005), the earthquake caused serious damage to oil storage tanks in Hokkaido, Japan. The earthquake provided a large number ofstrong-motion records, from which we can study how strong
ground motions were generated from the great subductionzone earthquake. In this study, a waveform inversion with
empirical Green’s functions was conducted to estimate
strong-motion generation areas (SMGAs) of the earthquakethat generated high-frequency ground motions up to 1 Hz.
The strong motions from the earthquake were observedin recently developed nationwide strong-motion networks inJapan including K-NET (Kinoshita, 1998) and KiK-net (Aoi
180

Strong-Motion Generation Areas of a Great Subduction-Zone Earthquake
Koketsu et al. (2004) conducted a joint inversion of strongmotion and geodetic data. These studies suggest that the rupture mainly propagated in the down-dip direction, with major
asperities away from the hypocenter.
Although these studies provided important constraints
on the rupture process at low frequencies, these studies effectively do not cover ground motions with frequencies higher
than 0.2 Hz. Although Yamanaka and Kikuchi (2003) andYagi (2004) used displacement waveforms with frequencies
lower than 0.5 or 1 Hz, their main focus is on frequencies
lower than 0.2 Hz as can be seen from figure 3 of Yamanaka
and Kikuchi (2003) or figure 8 of Yagi (2004). The restrictions on the frequency range in these studies arise from thedifficulty in computing realistic theoretical Green’s functions
at higher frequencies due to insufficient knowledge on subsurface structures in and around the source region.
In this study, we try to estimate strong-motion generation areas that produced ground motions with frequencies
from 0.1 to 1 Hz. Ground motions in this frequency range
are of engineering interest because of the following two reasons: First, many structures including tall buildings and longbridges have their natural frequencies in the same frequency
range. Second, sedimentary basins in Japan, which underlie
most of the major cities in Japan, tend to amplify groundmotions in this frequency range. It is therefore quite important to investigate how ground motions up to 1 Hz are generated from a great subduction-zone earthquake. This kind ofinvestigation will play a vital role in the effort to predict
strong ground motions from future huge subduction earthquakes including the scenario Tokai, Tonankai, and Nankai
earthquakes and to evaluate performance of structures subject to predicted ground motions.

181

Source processes for shallow crustal earthquakes havebeen successfully modeled in a deterministic approach fromthe waveform inversions using strong-motion data in the frequency range from 0.1 to 1 Hz (e.g., Miyakoshi et al., 2000).
Our intention is to reconstruct a fault rupture process effective in the same frequency range for the earthquake occurringin a subduction zone. We, therefore, use aftershock records
as empirical Green’s functions (Hartzell, 1978; Irikura, 1983)
instead of theoretical Green’s functions. In the empiricalGreen’s function method, it is necessary to use a small eventthat shares the path and the site effects with the large event.
Because the source region is large for the 2003 Tokachi-oki
earthquake, it is not appropriate to represent contributions
from the entire rupture with only one small event. A large
number of aftershocks following the 2003 Tokachi-oki earthquake were recorded by K-NET and KiK-net. We selected
several aftershocks to be used simultaneously in the inversion. Kamae and Kawabe (2004) also applied the empiricalGreen’s function method to the earthquake and developed asource model composed of three asperities with forward
modeling using records from a single aftershock. We selected
three aftershocks to be used in the analysis to represent theGreen’s functions from different parts of the mainshock rupture area more adequately.
In this study, we focus on the following two issues: The
first purpose is to confirm that the ground motions in thefrequency range up to 1 Hz from a great subduction-zone
earthquake can be modeled in a deterministic manner, once
realistic Green’s functions become available. The second
purpose is to compare the slip model obtained here withmodels so far obtained from ground-motion data at lower
frequencies. It is very important to investigate whether or

Figure 1. Observed velocity waveforms for the mainshock on the ground surface at TKCH07 (gray traces) and the calculated velocity
waveforms for Fourier amplitude of the mainshock and Fourier phase of aftershock 1 (black traces). Both traces are band-pass filtered
between 0.1 and 1.0 Hz. The similarity of the traces indicates that the mainshock and aftershock 1 share the path and site effects forthe site.

182

A. Nozu and K. Irikura

not SMGAs coincide with asperities that produce lowerfrequency ground motions. Such an investigation is essential in constructing reliable source models of future great
subduction-zone earthquakes for strong-motion prediction.

listed in Table 1; epicenters of the aftershocks are shown in
Fig. 2). Both traces are band-pass filtered between 0.1 and1 Hz. From Figure 1, we can recognize the similarity of the
traces.
The two traces (gray and black) in Figure 1 have thesame amplitude spectra as the mainshock but different phase
spectra: one is the mainshock, and the other is aftershock 1.
The group delay time (e.g., Sawada et al., 1998; Beauval
et al., 2003) is defined as the derivative of the phase spectrawith respect to angular frequency. The similarity of the two
traces means coincidence of the group delay time between
the mainshock and aftershock 1. In general, the group delay
time of a ground motion is the sum of the source, path, andsite effects as pointed out by Sawada et al. (1998) and Beauval et al. (2003). The source effect, however, is relatively
small as can be seen in the following example. Let us consider a simple source time function of an asperity represented
by a boxcar function with a duration of 2b (Fig. 3). Then the
Fourier transform of the function can be expressed as2=ω sinωb, which is a real number. Thus, both the Fourierphase and the group delay time of the function are zero for all
ω. Thus, the source process of the mainshock asperity does

Selection of Empirical Green’s Functions
In the empirical Green’s function method, it is important
to use appropriate Green’s functions; that is, it is necessary tofind a small event that shares the path and the site effects withthe large event. From the large number of aftershocks that
followed the 2003 Tokachi-oki earthquake, we selected appropriate aftershocks to be used in the analysis, referring to
the similarity of the Fourier phase or the “group delay time”
between the mainshock and the aftershock records asfollows.
In Figure 1, the gray traces indicate the observed velocity waveforms for the mainshock on the ground surface at
the KiK-net station TKCH07 (locations of the stations areshown in Fig. 2). On the other hand, the black traces indicate
calculated velocity waveforms for the Fourier amplitude of
the mainshock and the Fourier phase of aftershock 1 (source
parameters for the aftershocks as well as the mainshock are

45

HOKKAIDO ISLAND
Latitude (Deg.)

44
ABSH12

ABSH07
ABSH13
ABSH1543

HDKH02

NMRH05

NMRH04

KSRH06 KSRH04
KSRH10

KSRH07

TKCH05
TKCH06

TKCH11
HDKH01

KSRH05
KSRH02

TKCH03

TKCH04KKWH08

NMRH02
KSRH01
KSRH03

TKCH01
TKCH02

TKCH10

ABSH14

KSRH09

Aft.2

TKCH07
HKD096

IBUH03

HDKH05

IBUH05

TKCH08

HKD098

HDKH06

Aft.3
HDKH07

42

HKD113

Aft.1

141

142

143

144

145

146

Longitude (Deg.)

Figure 2.

A map showing the JMA hypocenters of the mainshock (open circle) and the aftershocks (solid squares), the fault plane
assumed for the inversion (large rectangle), the center of the rupture front assumed for the inversion (asterisk), and the K-NET and
KiK-net stations used for the inversion (solid triangles).

Strong-Motion Generation Areas of a Great Subduction-Zone Earthquake

183

Table 1
Parameters for the Mainshock and Aftershocks

Mainshock
Aftershock 1
Aftershock 2
Aftershock 3

Date (mm/dd/yyyy)

Time* JST

Latitude* (deg)Longitude* (deg)Depth* (km)MJ *M0 (Nm)09/26/2003
09/26/2003
09/27/2003
09/28/2003

4:50:07.64
7:20:22.49
17:06:21.56
9:23:11.30

41.78042.05742.73342.265144.078
143.734
144.346
143.324

42.0
40.9
59.2
43.4

8.0
5.2
5.2
5.0

8:21E20†
1:50E17‡
1:02E17†
3:95E16†(Strike, Dip, Rake) (deg)(249, 15, 127)†(216, 64, 87)§
(233, 26, 111)†(239, 20, 112)†*

After JMA.
After F-net (available at www.fnet.bosai.go.jp).
‡
Estimated from Fourier spectral ratios of aftershocks 1 and 2.
§
After Hi-net (available at www.hinet.bosai.go.jp).
†

not have significant effects on the group delay time of theobserved ground motion, although its effects on the Fourier
amplitude are significant. Thus, the similarity of the groupdelay time in Figure 1 suggests that the mainshock and aftershock 1 share the path and the site effects for the specific site.
More strictly, aftershock 1 shares the path and the site effects
with a part of the mainshock rupture that contributed most to
the ground motions at TKCH07. The similarity of the groupdelay time between the mainshock and aftershock 1 can beseen at many other stations. Because of such characteristics,
aftershock 1 is preferable as a small event to be used in theempirical Green’s function method. Aftershock 1 was also
used by Kamae and Kawabe (2004).
There are stations, however, where the mainshock groupdelay time is not similar to that of aftershock 1. The KiK-net
station KSRH09 (see Fig. 2), which is located east of
TKCH07, is such an example (Fig. 4). At KSRH09, thegroup delay time of the north–south component of the mainshock is different from that of aftershock 1. For this site, we
found that aftershock 2 (see Table 1 and Fig. 2) gives more
preferable results for both the east–west and the north–south
components (Fig. 5). In addition, we found a general tendency that at stations in the Kushiro region, that is, the stations with abbreviations beginning with KSR (see Fig. 2),
aftershock 2 gives more preferable results than aftershock 1.
At several stations in the west, neither aftershock 1 nor
aftershock 2 gives good results (Fig. 6), and after some trials
we found that aftershock 3 (see Table 1 and Fig. 2) gives the

Figure 3. A simple source time function of an asperity used toevaluate the effects of its rupture on the group delay time.

best results for those stations (Fig. 7). Based on such investigations, we decided to use aftershocks 1, 2, and 3 at thesame time in the ensuing analysis as illustrated in the next
section.
The same procedure to find appropriate aftershocks foruse in the empirical Green’s function method was also applied to the 2004 mid-Niigata Prefecture, Japan, earthquake
and contributed to the modeling of the near-source ground
velocities in the frequency range from 0.2 to 1 Hz
(Nozu, 2005).

Inversion Method
The stations used for the analysis are shown in Figure 2.
We mainly used the KiK-net borehole data (Aoi et al., 2000)
so that we can avoid significant effects of soil nonlinearity
during the mainshock. There were 39 KiK-net stations thatrecorded all aftershocks 1, 2, and 3 and the mainshock. Fourof the 39 stations, namely, HDKH04, ABSH06, ABSH11,
and SRCH09, were excluded from the analysis because someof the aftershock records did not show a satisfactory signalto-noise ratio. To improve azimuthal coverage in the nearsource region, we added 3 K-NET stations (Kinoshita,1998), namely, HKD096, HKD098, and HKD113, wheredata were obtained at the ground surface. The effect of soil
nonlinearity might be small, however, at these stations for thefrequency range of concern (0.1–1 Hz) because the thicknesses of the layers there with shear-wave velocity less than
300 m=sec are only 3, 1, and 1 m, respectively. As a result,the 38 stations shown in Figure 2 were used for the inversion.
The observed accelerograms for both the mainshock and theaftershocks were band-pass filtered between 0.1 and 1 Hz
and integrated in the frequency domain to obtain velocity
waveforms.
Radial components were used in the waveform inversion
considering the intensity of the radiation coefficient against
the diversity of the rake angle as follows. According to the
waveform inversion by Honda et al. (2004), the rake angle on
the mainshock fault plane is not necessarily uniform. In theempirical Green’s function method, however, the ground motions are synthesized assuming the same rake angle for themainshock and the aftershocks. If the radiation coefficient is

184

A. Nozu and K. Irikura

Figure 4. Observed velocity waveforms for the mainshock on the ground surface at KSRH09 (gray traces) and the calculated velocity
waveforms for Fourier amplitude of the mainshock and Fourier phase of aftershock 1 (black traces). Both traces are band-pass filtered
between 0.1 and 1.0 Hz. The discrepancy of the traces indicates that the mainshock and aftershock 1 do not share the path and site effects
for the site.
very sensitive to a small change in the rake angle, the waveform modeling with the empirical Green’s functions might
not be reasonable. We, therefore, investigated the sensitivityof the radiation coefficients for SV and SH waves to a smallchange in the rake angle for the mainshock. Figure 8 shows
the sensitivity of the radiation coefficients for SV and SH
waves to a small change in the rake angle defined as

Velocity (cm/s)

Δ

(1)

80

amplitude=mainshock, phase=mainshock
amplitude=mainshock, phase=aftershock2

KSRH09(EW)

40
0
−40
−80
0

Velocity (cm/s)

jR110  R90 j2;
jR110 jjR90 j

where R110 is the radiation coefficient for a point source at
the JMA hypocenter with a rake angle of 110° and R90 is thatwith a rake angle of 90°. The strike and dip angles were set tobe 246° and 18°, respectively (Honda et al., 2004). The radiation coefficients were calculated by applying ray tracing
assuming a layered subsurface structure used by Honda et al.
(2004). The results indicate that the radiation coefficients for
SV waves are much less sensitive to a small change in therake angle than those for SH waves in the areas wherestrong-motion sites are located. Thus, we decided that radial
components should be used in the waveform inversion. The

20

40

80

80

100

120

amplitude=mainshock, phase=mainshock
amplitude=mainshock, phase=aftershock2

KSRH09(NS)40

60Time (s)

0
−40
−80
0

20

40

60
Time (s)

80

100

120

Figure 5. Observed velocity waveforms for the mainshock on the ground surface at KSRH09 (gray traces) and the calculated velocity
waveforms for Fourier amplitude of the mainshock and Fourier phase of aftershock 2 (black traces). Both traces are band-pass filtered
between 0.1 and 1.0 Hz. The similarity of the traces indicates that the mainshock and aftershock 2 share the path and site effects forthe site.

Velocity (cm/s)

Strong-Motion Generation Areas of a Great Subduction-Zone Earthquake
20

0
−10
−20

Velocity (cm/s)

20

20

40

60
Time (s)

80

100

120

amplitude=mainshock, phase=mainshock
amplitude=mainshock, phase=aftershock1

HDKH01(NS)

10
0
−10
−20
0

Velocity (cm/s)

amplitude=mainshock, phase=mainshock
amplitude=mainshock, phase=aftershock1

HDKH01(EW)

10

0

20

20

40

60
Time (s)

80

100

120

amplitude=mainshock, phase=mainshock
amplitude=mainshock, phase=aftershock2

HDKH01(EW)

10
0
−10
−20
0

Velocity (cm/s)

185

20

20

40

80

100

120

amplitude=mainshock, phase=mainshock
amplitude=mainshock, phase=aftershock2

HDKH01(NS)1060Time (s)

0
−10
−20
0

20

40

60
Time (s)

80

100

120

Figure 6. Observed velocity waveforms for the mainshock on the ground surface at HDKH01 (gray traces) and the calculated velocity
waveforms for Fourier amplitude of the mainshock and Fourier phase of aftershock 1 (upper two black traces) or aftershock 2 (lower two
black traces). Both traces are band-pass filtered between 0.1 and 1.0 Hz. The discrepancy of the traces indicates that neither aftershock 1 nor
aftershock 2 shares the path and site effects with the mainshock for the site.
application of the present slip model to the transverse components will be discussed later.
The conventional least-squares linear waveform inversion (Hartzell and Heaton, 1983) was adopted. A fault planewith a dimension of 120 by 120 km was assumed (Fig. 2),
whose strike and dip angles were set to be 246° and 18°, respectively, to be consistent with the upper boundary of the
subducting Pacific plate (Honda et al., 2004). The fault
plane, that is, the rupture area, was divided into 30 by 30 elements. The rupture starting point is taken at a point (the asterisk in Fig. 2; 42.039° N, 143.925° E, 30-km depth) that is
approximately 30-km northwest of the JMA hypocenter. We
did not use the JMA hypocenter as the rupture starting point
because, as long as we use the JMA hypocenter, we could not
get good results. This is probably because the rupture front isassumed to propagate radially from the effective rupture

starting point in our analysis. The rupture is assumed to start
4 sec later than the JMA hypocentral time and to propagate
radially at a constant velocity of 2:6 km=sec. This value of
rupture velocity was selected based on the goodness of fit.
The relation between the rupture velocity and the goodnessof fit will be discussed later. Each fault element is allowed
to slip 12 times in 6 sec after passage of the rupture front at
0.5-sec time intervals. The moment release of each slip relative to the moment of the selected aftershock was determined through the inversion.
In terms of aftershocks, three cases of the inversions are
conducted (Fig. 9). In case 1, the contribution from the entire
rupture area is computed using records from aftershock 1. Incase 2, the contribution from the western part of the fault is
computed using records from aftershock 1, while that fromthe eastern part of the rupture area is computed using records

186
Velocity (cm/s)

A. Nozu and K. Irikura
20

0
−10
−20
0

Velocity (cm/s)

amplitude=mainshock, phase=mainshock
amplitude=mainshock, phase=aftershock3

HDKH01(EW)

10

20

40

20

80

100

120

amplitude=mainshock, phase=mainshock
amplitude=mainshock, phase=aftershock3

HDKH01(NS)1060Time (s)

0
−10
−20
0

20

40

60
Time (s)

80

100

120

Figure 7. Observed velocity waveforms for the mainshock on the ground surface at HDKH01 (gray traces) and the calculated velocity
waveforms for Fourier amplitude of the mainshock and Fourier phase of aftershock 3 (black traces). Both traces are band-pass filtered
between 0.1 and 1.0 Hz. The similarity of the traces indicates that the mainshock and aftershock 3 share the path and site effects forthe site.

Figure 8. Sensitivity of radiation coefficients for SV (left-hand panel) and SH (right-hand panel) waves to a small change in the rake angle
as defined in equation (1). If the sensitivity is larger than one, than it is set to be one in the figure. The radiation coefficients for SV waves aremuch less sensitive to a small change in the rake angle than those for SH waves in the areas where strong-motion sites are located.

Strong-Motion Generation Areas of a Great Subduction-Zone Earthquake

187

Figure 9. Three cases of inversions. In case 1, the contribution from the entire fault is computed using records from aftershock 1. Incase 2, the contribution from the western part of the fault is computed using records from aftershock 1, while that from the eastern part of thefault is computed using records from aftershock 2. In case 3, the contribution from the southwestern part of the fault is computed using
records from aftershock 1, the contribution from the eastern part of the fault is computed using records from aftershock 2, and the contribution
from the northwestern part of the fault is computed using records from aftershock 3.

from aftershock 2. In case 3, the contribution from the southwestern part of the rupture area is computed using records
from aftershock 1, that from the eastern part of the rupturearea is computed using records from aftershock 2, and thatfrom the northwestern part of the rupture area is computed
using records from aftershock 3.
Conventional corrections for the geometrical spreading
and time shifts (Irikura, 1983) were applied to the empiricalGreen’s functions to represent arrivals from each fault element. The aftershock locations determined by the JMA wereused (Table 1) in conjunction with the shear-wave velocity inthe source region of 3:8 km=sec (Ide and Takeo, 1996). Ab-

solute time information for both the mainshock and the aftershock records was used.
Nonnegative least-square solutions were obtained using
the algorithm of Lawson and Hanson (1974). In the analysis,
constraints were imposed to minimize the second order derivative of the slip on the rupture area. At boundaries of regions allotted to different aftershocks, the absolute momentfrom each element is assumed to be continuous, that is, therelative moment with respect to the aftershock is not necessarily continuous at the boundaries. The moments of aftershocks 2 and 3 are determined by the F-net (Fukuyama
et al, 1996) as listed in Table 1. The moment of aftershock

Figure 10. Slip models for the 2003 Tokachi-oki, Japan, earthquake obtained from three cases of waveform inversion. (a) The modelobtained from case 1 using only one aftershock; (b) the model obtained from case 2 using two aftershocks; and (c) the model obtained fromcase 3 using three aftershocks, which is the final inversion. The star indicates the rupture starting point.

188

A. Nozu and K. Irikura

1 was estimated from the averaged spectral ratio of aftershocks 1 and 2 at low frequencies at K-NET stations. In
Table 1, the mechanisms of the centroid moment tensor
solutions by the F-net for aftershocks 2 and 3 are listed.
For aftershock 1, the mechanism estimated from the polarity
distribution of the P wave by the Hi-net (Obara, 2002)
is shown.

Results
Inversion with One Aftershock

Velocity (cm/s)

Figure 10a shows the slip model obtained from the inversion case 1 using only one aftershock. The inversion
determines the distribution of the moment release in the rupture area relative to the aftershock moment. Then, using theaftershock moment listed in Table 1, absolute slip distribution in the rupture area is estimated as shown in Figure 10a,which is converted from the moment release distribution assuming a shear-wave velocity of 3:8 km=sec and a density of
3:0 g=cm3 (Ide and Takeo, 1996). The obtained slip model isquite simple. It has just one SMGA located close to the epicenter of aftershock 1. The synthetic velocities (0.1–1 Hz)
obtained from the model agree well with the observed one
at more than a few stations including TKCH07 (Fig. 11).
At some stations including KSRH03, however, the agreement is not satisfactory (Fig. 11). We, therefore, try to improve the result by including the second aftershock in theanalysis.

40

Inversion with Two Aftershocks
Figure 10b shows the slip model obtained from theinversion case 2 using two aftershocks. The obtained slip
model has two SMGAs. One is close to the epicenter of aftershock 1, which was also seen in the previous slip model, andthe other is located close to the northeastern corner of the
rupture area, closer to the epicenter of aftershock 2. The synthetic velocities (0.1–1 Hz) obtained from the model at threestations are compared with the observed ones in Figure 12.
The agreement between the synthetic and the observed velocities is again satisfactory for TKCH07. In addition, at
KSRH03, the agreement has been much improved.
The fitting between the synthetic and the observed velocities is evaluated by a residual, defined asv
u
N site R t2
u 1 X
t  Oi t2 dt
t
t1 S
Rit2
;
ERR 
2
N site i1t1 Oi t dt(2)

where Si t is the synthetic velocity at the ith site, Oi t is theobserved velocity at the ith site, N site is the number of sites
used for the inversion and t1 ; t2  is the time section used forthe inversion. As a result, the residual is 10% smaller forcase 2 than for case 1 (Table 2).
At some stations in the west including HDKH01 (seeFig. 2), however, the agreement is not satisfactory (Fig. 12).
We, therefore, try to further improve the result by including
the third aftershock in the analysis.

TKCH07

20
0
−20
−40
0

20

40

60

80

100

120

140

80

100

120

140

Velocity (cm/s)

Time (s)20

KSRH03

10
0
−10
−20
0

20

40

60
Time (s)Figure 11.

Comparison of the synthetic radial velocities (0.1–1 Hz) from the inversion case 1 (black traces) and the observed velocities(gray traces) at TKCH07 and KSRH03. Hatched portions are used for the inversion.

Velocity (cm/s)

Strong-Motion Generation Areas of a Great Subduction-Zone Earthquake
40

189

TKCH07

20
0
−20
−40
0

20

40

60

80

100

120

140

80

100

120

140

80

100

120

140

Velocity (cm/s)

Time (s)20

KSRH03

10
0
−10
−20
0

20

40

60

Velocity (cm/s)

Time (s)20

HDKH01

10
0
−10
−20
0

20

40

60
Time (s)Figure 12. Comparison of the synthetic radial velocities (0.1–1 Hz) from the inversion case 2 (black traces) and the observed velocies(gray traces) at TKCH07, KSRH03, and HDKH01. Hatched portions are used for the inversion.
Inversion with Three Aftershocks
Figures 10c and 13 show the slip model obtained fromthe final inversion (case 3) using three aftershocks. The obtained slip model has three SMGAs. The first one is close to
the hypocenter of aftershock 1, the second one is close to thatof aftershock 2, and the third one is close to that of aftershock 3 (Fig. 13). These SMGAs are hereafter referred toas SMGA 1, SMGA 2, and SMGA 3. The total moment of
the model is 2:06E21 Nm (Mw 8:1), which includes contributions from the SMGAs and the background area. The
synthetic velocities (0.1–1 Hz) obtained from the model at
all stations are compared with the observed ones in Figure 14.
The agreement is satisfactory at almost all the stations including TKCH07, KSRH03, and HDKH01. The residual defined by equation (2) is 20% smaller for case 3 than for case 1(Table 2). Thus, we prefer the slip model from case 3. The
goodness of fit is mainly due to the use of empirical Green’s
functions, which allows us to represent complicated path andsite effects in the frequency range from 0.1 to 1 Hz.

The reason why case 1 and case 2 cannot estimate someof the asperities can be explained as follows. Because thesource region of the 2003 Tokachi-oki earthquake is so large,
it is not appropriate to represent contributions from the entire
rupture area with only one aftershock. For example, recordings from aftershock 1 are only appropriate for representing
arrivals from the southwestern part of the fault; they do notshare the path and the site effects with arrivals from otherparts of the fault. Thus, as long as we use records from only
aftershock 1, we cannot achieve good fits for parts of theTable 2
ERR Values Defined in Equation (2)for Three Inversion Cases
Case

ERR

1
2
3

0.617
0.555
0.491

190

Figure 13.

A. Nozu and K. Irikura

Slip model for the 2003 Tokachi-oki, Japan, earthquake obtained from the final inversion (case 3) using three aftershocks.

mainshock records that are affected by SMGA 2 or SMGA 3.
This is why case 1 cannot estimate SMGA 2 and SMGA 3 andcase 2 cannot estimate SMGA 3.
Figure 15 shows the slip for each time window of 4 sec
during the 2003 Tokachi-oki earthquake obtained by the final
inversion. It can be observed that the rupture of SMGA 1 took
place from 4 to 12 sec, while the rupture of SMGA 2 took
place from 32 to 40 sec. The rupture of SMGA 3 took placefrom 20 to 28 sec, although the slip velocity is rather small
compared to SMGA 1 and SMGA 2. Figure 16 shows representative slip velocity time functions for each SMGA. The
time functions correspond to the element with the largest
final slip in each SMGA. The slip velocity should be regarded
as the averaged value over each time window used for the
waveform inversion (0.5 sec). In SMGA 1, where the largest
slip occurred, the peak slip velocity averaged over the time
window of 0.5 sec is estimated to be at least as highas 3 m=sec.
Figure 17 shows the relation between the rupture velocity and ERR in equation (2). Based on the relation, we preferred the value of 2:6 km=sec as the rupture velocity.

Discussion
Locations of SMGAs and Selected Aftershocks
In the preceding analysis, we selected three aftershocks
for use in the analysis, referring to the similarity of the groupdelay time between the mainshock and the aftershock records. The inversion result indicates that the mainshock is
composed of three SMGAs, each of which is close to the location of one of the aftershocks (Fig. 13). From this configuration of the SMGAs and the selected aftershocks, we can
explain the reason for the similarity of the group delay time
between the mainshock and the aftershock records (Figs. 1,
5, and 7) as follows:
1. Ground motions at a site from the mainshock are roughly
expressed as a sum of contributions from the three
SMGAs. Some sites are affected dominantly by one of
the three SMGAs and less by the other two, mostly depending on distance between the site and the SMGA.
2. Then, the observed records of an aftershock close to the
SMGA are appropriate for the empirical Green’s functions, because the path and site effects are shared between

Strong-Motion Generation Areas of a Great Subduction-Zone Earthquake

Figure 14.

191

Comparison of the synthetic radial velocities (0.1–1 Hz) from the final inversion (case 3) (red traces) and the observed
velocities (black traces) at all stations used for the inversion. The stations are aligned largely from west to east. Hatched portions are used
for the inversion.
(Continued)

192

A. Nozu and K. Irikura

Figure 14.

Continued.

Strong-Motion Generation Areas of a Great Subduction-Zone Earthquake

Figure 14.
the mainshock and the aftershock, resulting in a similarityof the group delay time between the mainshock and theaftershock for the specific site.

Comparison with Low-Frequency Slip
Models So Far Obtained
As mentioned before, it is very important to investigate whether or not strong ground motions from a great
subduction-zone earthquake are generated from the asperities defined as large slip areas from the waveform inversionof lower-frequency ground motions. Generally speaking, locations of the SMGAs identified in this study agree well withthe asperities identified from lower-frequency ground motions in the previous studies.

193

Continued.

The location of SMGA 1 in this study (Fig. 13) is almost
the same as the location of the major asperity identified fromthe joint inversion of strong-motion and geodetic data by Koketsu et al. (2004). It is also almost the same as the region ofmajor uplift identified from first-peak tsunami travel times(Hirata et al, 2004; region A of their fig. 3). The major asperity by Yamanaka and Kikuchi (2003), the asperity B by
Yagi (2004) and the southern half of the asperity (b) by
Honda et al. (2004) are close to our SMGA 1. SMGA 2 inthis study (Fig. 13) is close to the asperity (c) by Honda
et al. (2004). Areas of relatively large slip can also be foundin the models of Yamanaka and Kikuchi (2003), Yagi (2004),
and Koketsu et al. (2004), close to our SMGA 2, although
they are not referred to as asperities in those studies. SMGA 3
in this study (Fig. 13) is close to the asperity C by Yagi

194

A. Nozu and K. Irikura

Figure 15.

Slip for each 4 sec of the 2003 Tokachi-oki earthquake obtained by the final inversion.

(2004) and the northern half of the asperity (b) by Honda
et al. (2004).
Thus, the SMGAs identified in this study agree well withasperities identified from inverted results so far obtainedusing lower-frequency ground motions plus geodetic data
and tsunamis. It implies that strong ground motions up to
1 Hz were generated from almost the same asperities producing lower-frequency ground motions. In terms of the complexity of slip models, although our analysis is focused onhigh frequencies, our slip model is at least as simple asconventional low-frequency slip models. Such results would
be useful in constructing source models of future great
subduction-zone earthquakes for strong-motion prediction.
The rupture velocity preferred in our analysis(2:6 km=sec) is smaller than those used in the conventional
low-frequency inversions, for example, 4:5 km=sec suggested by Yagi (2004). This is probably because our analysis
is focused on the main rupture that generated high frequencies following the rupture of SMGA 1, and thus, our value of
rupture velocity is mainly controlled by the time lag between
the rupture of SMGA 1 and that of SMGA 2. Kamae and Ka-

wabe (2004), who are also focussed on high frequencies,
used a relatively small rupture velocity (2:8 km=sec), whichis comparable to ours.
Problems in Transverse Components
Finally, the limitation of the present model associatedwith the applicability to the transverse component is discussed. Figure 18 shows comparisons of observed (graytraces) and synthetic (black traces) transverse ground velocities, band-pass filtered between 0.1 and 1 Hz, at some stations. The synthetic velocities were computed from the final
slip model (Fig. 13). At TKCH05, TKCH08, and KSRH09,
transverse synthetic ground velocities agree well with the observed ones. At TKCH03, however, transverse ground velocity is underestimated. As mentioned earlier, the radiation
coeffieients for SH waves are much more sensitive to a smallchange in the rake angle than those for SV waves in the areas
where strong-motion sites are located (Fig. 8). Such a sensitivity may have caused the discrepancy found in the transverse components.

Strong-Motion Generation Areas of a Great Subduction-Zone Earthquake

195

0.70

4

0.65

3

0.60

2

ERR

Slip velocity (m/s)

SMGA 1

1

0.55
0.50

0
0

1

2

3

4

5

6

Time after passage of rupture front (s)40.40Slip velocity (m/s)

SMGA 2

2.0

3

2.5 3.0 3.5 4.0 4.5
Rupture Velocity (km/s)

2

Figure 17. Relation between the rupture velocity and ERR inequation (2). ERR is minimized when we use a rupture velocityof 2:6 km=sec.

1

0
0

1

2

3

4

5

6

Time after passage of rupture front (s)4SMGA 3
Slip velocity (m/s)

0.45

3

2

1

0
0

1

2

3

4

5

6

Time after passage of rupture front (s)Figure 16.

Representative slip velocity time functions for each

SMGA. The time functions correspond to the element with the largest final slip in each SMGA. The slip velocity should be regarded as

the averaged value over the time window used for the waveform
inversion (0.5 sec).

Conclusions
A waveform inversion was conducted for the 2003
Tokachi-oki, Japan, earthquake (MJ 8.0) using empiricalGreen’s functions to estimate SMGAs of the earthquake thatgenerated high-frequency ground motions with frequencies
up to 1 Hz. We selected three aftershocks for use in theanalysis, referring to the similarity of the group delay time
between the mainshock and the aftershock records. Threecases of inversions are conducted, namely, case 1 in which
only one aftershock was used, case 2 in which two aftershocks were used, and case 3 in which three aftershocks wereused. Case 3 yielded the best results. The resultant slip modelhas three SMGAs, each of which is close to the location of
one of the aftershocks used as an empirical Green’s function
event. We found that careful selection of aftershocks is very
effective for estimating ground motions from the mainshock,taking into account the similarity of the group delay time
between the mainshock and the aftershock records. We confirmed that the ground motions with frequencies up to 1 Hz
from a great subduction-zone earthquake can be modeled in a
deterministic manner, once realistic Green’s functions become available. Generally speaking, locations of the SMGAs
identified in this study agree well with asperities identified
from inverted results using lower-frequency (lower than
0.2 Hz) ground motions plus geodetic data and tsunamis.
It implies that strong ground motions up to 1 Hz were generated from the asperities producing lower-frequency groundmotions. In terms of the complexity of slip models, although
our analysis is focused on high frequencies, our slip model is
at least as simple as conventional low-frequency slip models.
Such results would be useful in constructing source models
of future great subduction-zone earthquakes for strongmotion prediction.

Acknowledgments
We would like to express our gratitude to the National Research Institute for Earth Science and Disaster Prevention (NIED), Japan, for providing important strong-motion data from K-NET and KiK-net and the moment
tensor catalog from F-net. We would also like to thank the Japan Meteorological Agency for providing source parameters for the 2003 Tokachi-oki

196

Velocity (cm/s)

A. Nozu and K. Irikura
20

TKCH03

10
0
−10
−20
0

20

40

60

80

100

120

140

80

100

120

140

80

100

120

140

80

100

120

140

Velocity (cm/s)

Time (s)20

TKCH05

10
0
−10
−20
0

20

40

60

Velocity (cm/s)

Time (s)20

TKCH08

10
0
−10
−20
0

20

40

60

Velocity (cm/s)

Time (s)KSRH09

40
0
−40
0

20

40

60
Time (s)Figure 18.

Synthetic transverse velocities (black traces) from the final slip-model (Fig. 13) are compared with the observed velocities(gray traces) at four stations. At TKCH03, the observed transverse velocity is characterized by a large-amplitude pulselike motion, which is
not captured by the synthetic velocity.

mainshock and aftershocks. The manuscript was significantly improved by
the comments from Hiroshi Kawase, Kyushu University, and two anonymous reviewers.

References
Aoi, S., K. Obara, S. Hori, K. Kasahara, and Y. Okada (2000). New strongmotion observation network: KiK-net, EOS 81, 329.
Beauval, C., P. Y. Bard, P. Maczo, and J. Kristek (2003). Quantification of
frequency-dependent lengthening of seismic ground motion duration
due to local geology: application to the Volvi Area (Greece), Bull.
Seismol. Soc. Am. 93, 371–385.

Fukuyama, E., M. Ishida, S. Hori, S. Sekiguchi, and S. Watada (1996).
Broadband seismic observation conducted under the FREESIA Project, Rep. Natl. Res. Inst. Earth Sci, Disaster Prev. 57, 23–31.
Hartzell, S. H. (1978). Earthquake aftershock as Green’s functions, Geophys.
Res. Lett. 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. Seismol. Soc. Am. 73,
1553–1583.
Hatayama, K., S. Zama, H. Nishi, M. Yamada, Y. Hirokawa, and R. Inoue(2004). Long-period strong ground motion and damage to oil storage
tanks due to the 2003 Tokachi-oki earthquake, J. Seism. Soc. Jpn. 57,
83–103 (in Japanese).
Hirata, K., Y. Tanioka, K. Satake, S. Yamaki, and E. L. Geist (2004). The
tsunami source area of the 2003 Tokachi-oki earthquake estimatedStrong-Motion Generation Areas of a Great Subduction-Zone Earthquake
from tsunami travel times and its relationship to the 1952 Tokachi-oki
earthquake, Earth Planet. Space 56, 367–372.
Honda, R., S. Aoi, N. Morikawa, H. Sekiguchi, K. Kunugi, and H. Fujiwara(2004). Ground motion and rupture process of the 2003 Tokachi-oki
earthquake obtained from strong motion data of the K-NET and
KiK-net, Earth Planet. Space 56, 317–322.
Ide, S., and M. Takeo (1996). The dynamic rupture process of the 1993
Kushiro-oki earthquake, J. Geophys. Res. 101, 5661–5675.
Irikura, K. (1983). Semi-empirical estimation of strong ground motions
during large earthquakes, Bull. Disaster Prev. Res. Inst. 32,
63–104.
Kamae, K., and H. Kawabe (2004). Source model composed of asperities forthe 2003 Tokachi-oki, Japan, earthquake (MJMA  8:0) estimated by
the empirical Green’s function method, Earth Planet. Space 56,
323–327.
Kinoshita, S. (1998). Kyoshin-net (K-NET), Seism. Res. Lett. 69, 309–332.
Koketsu, K., K. Hatayama, T. Furumura, Y. Ikegami, and S. Akiyama(2005). Damaging long-period ground motions from the 2003
Mw 8.3 Tokachi-oki, Japan, earthquake, Seism. Res. Lett. 76, 67–73.
Koketsu, K., K. Hikima, S. Miyazaki, and S. Ide (2004). Joint inversion ofstrong motion and geodetic data for the source rupture process of the2003 Tokachi-oki, Hokkaido, earthquake, Earth Planet. Space 56,
329–334.
Lawson, C. L., and R. J. Hanson (1974). Solving Least Squares Problems,
Prentice Hall, Englewood Cliffs, New Jersey.
Miyakoshi, K., T. Kagawa, H. Sekiguchi, T. Iwata, and K. Irikura (2000).
Source characterization of inland earthquakes in Japan using source
inversion results, in Proc. 12th World Conf. Earthquake Engineering, Aukland, New Zealand, 30 January–4 February: available on
CD-ROM.
Nozu, A. (2005). Variable-slip rupture model for the 2004 Mid Niigata Prefecture (Niigata-ken Chuetsu) earthquake—waveform inversion with

197

empirical Green’s functions, Zisin 58, 329–343 (in Japanese with English abstract).
Obara, K. (2002). Hi-net: high sensitivity seismograph network, Japan, Lect.
Note Earth Sci. 98, 79–87.
Sawada, S., H. Morikawa, K. Toki, and K. Yokoyama (1998). Identificationof path and local site effects on phase spectrum of seismic motion, in
Proc. 10th Japan Conf. Earthquake Engineering, Yokohama, Japan,
25–27 November, 915–920.
Tsuchida, H., E. Kurata, and K. Sudo (1969). Technical note no. 80, Strongmotion earthquake records on the 1968 Tokachi-oki earthquake and its
aftershocks, Port and Harbour Research Institute (in Japanese).
Yagi, Y. (2004). Source rupture process of the 2003 Tokachi-oki earthquake
determined by joint inversion of teleseismic body wave and strong
ground motion data, Earth Planet. Space 56, 311–316.
Yamanaka, Y., and M. Kikuchi (2003). Source process of the recurrent
Tokachi-oki earthquake on September 26, 2003, inferred from teleseismic body waves, Earth Planet. Space 55, e21–e24.
Geotechnical and Structural Engineering Department
Port and Airport Research Institute3-1-1, Nagase
Yokosuka 239-0826, Japan(A.N.)
Disaster Prevention Research Center
Aichi Institute of Technology
1247 Yachigusa
Yagusa, Toyota 470-0392, Japan(K.I.)
Manuscript received 28 August 2006