SOURCE RUPTURE PROCESS OF THE MENTAWAI,INDONESIA EARTHQUAKE OF OCTOBER 25th, 2010,
DETERMINED BY JOINT INVERSION OFTELESEISMIC BODY WAVE AND NEAR-SOURCESTRONG MOTION DATA
Iman Fatchurochman
MEE10502

Supervisor:Yuji YAGIABSTRACT
We investigated the source process of the 2010 Mentawai earthquake by joint inversion method usingnear source and teleseismic body-wave data. To perform a stable inversion, we applied smoothing
constraints and determined their relative weights on the observed data using ABIC criterion. Theteleseismic waveforms were windowed for 150 sec, band-passed between 0.001 and 1.0 Hz, and thenintegrated into displacement. The strong motion data were windowed for 250 sec, band-passed
between 0.005 and 0.5 Hz, and then integrated into velocity. The re-sampling rate for both data set tobe 0.25 sec. We estimated the fault area to be 190 × 70 km2. The main source parameters are as
follows: (strike, dip, rake) = (324⁰, 10⁰, 94.6⁰); the depth of hypocenter 13.5 km; the seismic moment
5.814×1020 Nm (Mw 7.8); source duration 102 sec; and the maximum slip amounts to 3.9 m.
The source rupture process obtained probably can be divided into 2 stages. At stage 1, therupture nucleated near the hypocenter and then propagated to the southwestward and broke the firstasperity centering at 14 km from the epicenter. At stage 2, the rupture propagated to the northwestward
and broke the second asperity which was centered about 78 km from the epicenter. Our total slip
distribution is well consistent with the result of tsunami waveform analysis by Satake et al. (2011).
This earthquake was categorized as a tsunami earthquake due to the long rupture duration and thegeneration of tsunami much larger than expected for its magnitude.
Keywords: 2010 Mentawai Earthquake, Source Process, Tsunami Earthquake, Joint Inversion.

1. INTRODUCTIONOn October 25th, 2010, a tsunami earthquake with a moment magnitude of 7.7 occurred closed toMentawai Island, Indonesia. The location of this earthquake is near the trench where theIndo-Australia Plate is subducted beneath the Eurasia Plate. Some researchers identify this earthquake
as a tsunami earthquake because the duration of the rupture was excessive and it generated a tsunamithat was much larger than expected for its seismic magnitude (Newman et al. 2011, Lay et al. 2011).
According to USGS data, the Mentawai earthquake has moment magnitude Mw 7.7 and surface
magnitude Ms 7.2. This disparity of Ms-Mw is also characteristic of tsunami earthquake (Kanamoriand Kikuchi, 1993). We can classify the earthquake that generates a tsunami as tsunamigenic, while
the tsunami earthquake is a special class of events that generate tsunamis much larger than expected
for their surface waveform magnitude (Kanamori, 1972). To understand the mechanism of tsunami

Indonesian Agency for Meteorology, Climatology and Geophysics (BMKG).
Professor, Graduate School of Life and Environmental Sciences, University of Tsukuba, Japan.

1earthquake, it is important to investigate the coseismic slip area of these great earthquakes using local
and global seismic data set.

2. DATAIn this study, 33 teleseismic stations (IRIS-DMC) and 12 components of five strong motion stations(BMKG) were used to infer the source process of the Mentawai earthquake 2010. The teleseismic
body wave waveforms were windowed for 150 sec with starting 10 sec before P-wave arrival, band
pass filtered between 0.001 and 1.0 Hz, and then integrated into displacement with sampling time 0.25
sec. The near source data were windowed for 250 sec with starting 5 sec before P-wave arrival andband pass filtered between 0.005 and 0.5 Hz, and then integrated into velocity with sampling time 0.25
sec.

b)

a)

Sumatra Island
Mentawai Island

Figure 1. Station distribution is shown as a map view. The star represents the epicenter of themain-shock, red tringles indicate the stations. a) teleseismic network, b) strong motion network.

3..THEORY AND METHODOLOGY
To obtain the source process of the Mentawai earthquake, we analyzed the event using the standard
fault parameterization. An inversion method was developed to infer earthquake rupture process fromfar-field and near source data set. The basic theory follows the standard waveform inversion scheme(Hartzell and Heaton, 1983) and numerical method developed by Yagi et al. (2004). The rupture
process as a spatiotemporal slip distribution is represented on a fault plane. The fault plane is divided
into M×N subfaults with length dx and width dy. Let xm= m dx and yn = m dy. Slip rate functions oneach subfault are described by a series of triangle functions with rise time τ and the fault slip vector
described by K basis slip vectors. Using this source model, the observed seismic waveform at thestation j is stated by:
𝑊𝑗𝑜𝑏𝑠 𝑡𝑖 =𝑋𝑚𝑛𝑙𝑘 𝑔𝑚𝑛𝑘𝑗 𝑡𝑖 − 𝑙 − 1 𝜏 − 𝑇𝑚𝑛 + 𝑒𝑗

(1)

𝑚𝑛𝑙𝑘𝑇𝑚𝑛 =2 +𝑦 2𝑥𝑚𝑛𝑉(2)where Xmnkl is the kth component of slip at the mnth subfault with the lth time step; gmnkj(t) is the green’sfunction at the mn-subfault; Tmn is the start time of the basis function at each subfault and 𝑉 is therupture velocity; ej is the Gaussian error with variance of j. We calculate Green’s function of2

subfaults assuming the point source at thecenter of the subfaults. The sampling time of
the Green’s function was set at 0.25 sec, samewith the observed waveform. The velocity
structures assumed for theoretical teleseismic
body wave and near source data calculation are
based on the crustal structure of the central
Sunda margin at the onset of obliquesubduction (Kopp et al. 2001) as shown intable 1.
Generally, an increase in the numberof model parameters may lead to give
instability in the solution even a small change
in the data resulting in a large change in thesolution (Yagi et al. 2004). Considering the instability in the solution, we applied smoothing
constraints to the slip distribution with respect to time and space to obtain more stable results. Thesmoothness constraint with respect to time is:Table 1. Velocity structure (Kopp et al. 2001).

0 = 𝑋𝑚𝑛𝑘 (𝑙−1) − 2𝑋𝑚𝑛𝑘𝑙 + 𝑋𝑚𝑛𝑘 (𝑙+1) + 𝑒𝑡,

2 ≤ l ≤ L-1(3)

where et is the Gaussian error. We can rewrite the equation (3) in the following simple vector form
0 = Tx +et,

(4)

where T is N1 x Na matrix (N1 = MN(L-2)K). The smoothness constraint on the spatial distribution of
total slip or smoothing constraint respect to space is introduced by a Laplacian finite difference
operator:
0=𝑋(𝑚 −1)𝑛𝑘𝑙 + 𝑋 𝑚 +1 𝑛𝑘𝑙 + 𝑋𝑚 (𝑛−1)𝑘𝑙 + 𝑋𝑚 𝑛+1 𝑘𝑙 − 4𝑋𝑚𝑛𝑘𝑙 + 𝑒𝑑𝑙(5)𝑋0𝑛𝑘𝑙 = 𝑋𝑚 0𝑘𝑙 = 𝑋(𝑀+1)𝑛𝑘𝑙 = 𝑋𝑚 (𝑁+1)𝑘𝑙 = 0(6)The equation (5) can be rewritten as follows:
0 = Dx +ed,

(7)

where D is N2 x Na matrix (N2 = MNK). To determine the model parameters that minimize the sum ofsquare residual S from the observed data with two constraints, given by
𝑆 𝒙, 𝜎𝒋 , 𝜎𝒕 , 𝜎𝒅 =𝑗1
1
12+ 2 𝑻𝒙 2 + 2 𝑫𝒙 22 𝒚𝑗 − 𝑨𝑗 𝒙𝜎𝑗𝜎𝑡𝜎𝑑(8)The values of t and d cannot be estimated directly, but j can be estimated by the quality of data. Todetermine t and d objectively, we adopted the minimum Akaike’s Bayesian information criterion(ABIC) (Akaike, 1980). The optimal ABIC for the present case is expressed using the following
equation:
𝐴𝐵𝐼𝐶 𝒙, 𝜎𝒕 , 𝜎𝒅 = 𝑵 log 𝑆 𝒙, 𝜎𝒋 , 𝜎𝒕 , 𝜎𝒅 − log

1𝜎𝒅𝟐𝑫𝒕𝑫+C

1𝜎𝒕2

1𝑻𝒕 𝑻 + 𝜎 𝟐 𝑫𝒕 𝑫 + log
𝒅1 𝑡
1𝒕j 𝜎 2 𝐴𝑗 𝐴𝑗 + 𝜎 2 𝑻 𝑻 +𝑗𝒕(9)N is the total number of the observation equations. We apply a grid-search method to obtain optimalvalues of t and d. We employed the non-negative least squares (NNLS) algorithm of Lawson and3

Hanson (1974) to give positivity constraints on the model parameter. The positivity constraint is
imposed not only because it is physically reasonable but also because negative slips lead to destructive
interference between subfaults, producing unstable solutions (Hartzell and Heaton, 1983).

4. RESULTS AND DISCUSSION
4.1. Finite Fault Model
Assuming that the faulting occurs on the single fault
plane, we constructed the fault plane according to theaftershock distribution in a week after the mainshock andadjusted it to extend near to the seafloor as well. We thendivided the fault plane into many subfaults and adopted
the geometry of the focal mechanism from global CMT
solution (strike, dip, rake) = (316⁰, 8⁰, 96⁰) with slight
modification to be consistent with the geometry of theplate boundary (Bird, 2003) as this result (strike, dip,rake) = (324⁰, 10⁰, 96⁰) with the epicenter determined byBMKG (latitude = -3.49⁰; longitude = 100.14⁰). Figure 2
shows the geometry of the fault plane where the redrectangle denotes the geometry of global CMT solution
result, while the blue rectangle denotes the geometry after
Figure 2. Geometry of the fault plane
has slightly been modified from global being modified.
To obtain stable and detailed rupture process, we
CMT solution result.
tried many lengths and widths of the fault plane. We
finally constructed the fault plane with the area of 190 km × 70 km, the dimension of subfault 10 km2
and the dip angle of 10⁰ to cover the slip distribution of this earthquake, after many tries and errors.
Assuming that the earthquake occurred in interplate and considering the structure model andthe dip angle, we adjusted the fault plane to extend near to the surface and we estimated the depth of
hypocenter to be 13.5 km. The slip rate function of the subfault is expanded into a series of 15 triangle
functions with a rise time of 2 sec. We also examined the optimum value for Vi within the range of 0.5
km/sec to 3.5 km/sec, and found that the maximum rupture velocity was 2.0 km/sec, which showed aminimum variance. The rigidity around the fault plane was assumed to be μ = 26.81 GPa from thestructure model (Table 1).
4.2. Source Rupture Process
The inversion results are shown in Figure 3. Figures 3 (a), (b), (c) and (d) show the fault type, thesource duration, the slip distribution and the comparison of observed with calculated waveforms ofthis earthquake, and Figure 4 shows the map view of slip distribution. Figure 3(a) illustrates the thrustfault type for this earthquake. The source time function (Figure 3b) indicates that a total seismic
moment M0 of 0.5814×1021 Nm (Mw 7.8) was released during a period of 102 sec. The total seismic
moment estimated by this study is in agreement with global CMT moment solution, 0.677×1021 Nm(Mw 7.8). This result shows the anomalously long rupture duration comparing with those ordinary
earthquakes with similar magnitude, such as the 1996 Peru earthquake (Mw 7.6 – 7.7), whose rupture
duration is around 45 to 50 sec (Jennifer et al. 1999). Due to the excessive long rupture duration thatreached 102 sec and generation of tsunami much larger than expected for its magnitude, the 2010
Mentawai earthquake is categorized as a tsunami earthquake. The final dislocation (Figure 3c)
describes that the slip distribution extends along the dip direction. The rupture propagated from thehypocenter to the shallower subfaults. In asperity 1, the maximum slip is around 3.9 m, while inasperity 2 it is around 2.5 m.

4Figure 3. Joint inversion result of the Mentawai earthquake. a)
Focal mechanism, b) source time function, c) slip distributionand d) comparison between observed (black) and calculated(red) waveforms

Figure 4. Distribution of coseismic
slip on the map and aftershocks
distribution in a week

This result is consistent with the maximum slip more than 3.0 m found by the result of tsunami
waveform analysis (Satake et al. 2011). Our inversion result shows that in the shalower part, the slip
ranges from 0.8 to 2.5 m, which is in the range of tsunami inversion result obtained by Satake et al.
(2011). According to this proof, it probably can explain why this earthquake generated tsunami.
However, our total slip distribution is different from the result of tele-seismic body wave analysis by
Newman et al. (2011). They found the rupture was centered in one asperity with the maximum slip
which was scaled from the original result to be 10 m, whereas our total slip distribution is wellconsistent with the result of tsunami waveform analysis by Satake et al. (2011).
Figures 3(d) shows the comparison between the observed (black curve) and the synthetic
waveforms (red curve). Generally the waveforms show a good agreement between them except thatthe MNSI station located around 400 km north of the hypocenter is not fitting well. The discrepancy
for this station may come from the assumed structure, which is not suitable for this station. According
to this result, in general the assumed structure model we used seems to be suitable for all stations.
Figure 4 shows the map view of coseismic slip and aftershock distribution in a week after themainshock. The aftershock data are obtained from BMKG while the focal mechanisms are taken fromglobal CMT solution results. The rupture propagates mainly along the dip direction and the rupturearea appears to extend parallel to the subduction zone between the Indo-Australia and Eurasia Plates.
The focal mechanisms of the aftershocks as shown in Figure 4 predominantly are normal faulting near
the trench. It probably describes the rupture of a locked subducted seamount on an otherwise
decoupled zone, resulting in extension of the outer‐rise causing the normal faulting aftershocks.
From the azimuthal view, the near source data only cover the epicenter from one side.
Although this data set is not quite good especially in azimuthal coverage context, our result seems tohave a good resolution. Probably this result may come from the advantage of the joint inversion ofteleseismic body wave and strong motion data with an optimized ABIC. Fukahata et al. (2003) showsthat optimized ABIC is particularly useful even with the insufficient observed data set.
The source rupture process obtained is probably divided into 2 stages. At stage 1, the rupture
nucleated near the hypocenter and then propagated to the southwest direction and broke the firstasperity centering at 14 km from the epicenter during 4 to 20 sec after the initial break. At this stage
the maximum slip amounts to 3.9 meter at 20 sec after the initial break with the seismic moment of
2.359 × 1020 Nm (Mw = 7.5) which is smaller than the total seismic moment. At stage 2, the rupture
propagated to the northwest direction and broke the second asperity which was centered about 78 kmfrom the epicenter. The seismic moment for this asperity is 2.368 × 1019 Nm (Mw = 6.9).

5

5. CONCLUSIONS
By Joint Inversion method using the near source strong motion data obtained by BMKG andteleseismic data collected by IRIS-DMC, we estimated a detailed source model and slip distribution of
the Mentawai, Indonesia earthquake of October 25th 2010. To perform a stable inversion, we appliedsmoothing constraints and determined their relative weights on the observed data using ABIC.
Assuming that the faulting occurs on the single fault plane, we constructed the fault plane
according to the aftershock distribution in a week after the mainshock and adjusted it to extend near tothe surface as well. We estimated the fault area to be 190 × 70 km2. The main source parameters are as
follows: (strike, dip, rake) = (324⁰, 10⁰, 94.6⁰) and the depth of hypocenter 13.5 km. We obtained total
seismic moment Mo = 0.5814x1021 Nm correspond to Mw = 7.8 with maximum slip amounting to 3.9m at 14 km southwest of the epicenter. Source duration of this earthquake is around 102 sec.
The rupture process obtained is probably divided into 2 stages: the rupture nucleated aroundthe hypocenter and propagated to the southwest and broke the first asperity centering at 14 km fromthe epicenter with maximum slip amounting to 3.9 m (stage I), then the rupture propagated to thenorthwest and the second asperity was broken, which was centered about 78 km from the epicenter(stage II). The rupture propagated mainly along the dip direction. We identify this earthquake as atsunami earthquake in that it has characteristic of excessive long rupture duration and generated atsunami that was much larger than expected for its magnitude. Our total slip distribution is wellconsistent with the result of tsunami waveform analysis by Satake et al. (2011).

6. RECOMMENDATION
The determination of source rupture process of tsunami earthquake is relatively difficult than ordinary
earthquake. Some researchers found the result of slip distribution obtained by seismic inversion
analysis was under estimation, comparing with the tsunami field observation. Therefore, we had better
evaluate the slip distribution obtained by seismic inversion analysis as input for tsunami simulation
model. We can explain the consistence of the seismic inversion analysis by comparing the synthetic
waveforms obtained by tsunami simulation with the observed tsunami waveforms.
AKNOWLEDGEMENT
I would like to express my sincere gratitude to my advisor Dr. B. Shibazaki for his continuous support,
valuable suggestion and guidance during my study.

REFERENCES
Akaike, H., 1980, pp. 143-166, University Press, Valencia, Spain.
Bird, P., 2003, Geochemistry Geophysics Geosystems, 4(3), 1027, doi:10.1029/2001GC000252.
Fukahata, Y., Y. Yagi., and M. Matsu'ura., 2003, Geophys. Res. Lett., 30, 10.1029/2002GL016293.
Hartzell, S. H. and T. H. Heaton., 1983, Bull. Seism. Soc. Am., 73, 1553-1583.
J. L. Swenson and S. L. Beck., 1999, Pure appl geophys Volume 154, Numbers 3-4, 731-751
Kanamori, H., 1972, Phys.Earth Planet. Interiors, 6, 346-359.
Kanamori, H., M. Kikuchi., 1993, Nature, 361.
Kopp, H., et al., 2001, Geophys. J. Int, 147, 449-474.
Lawson, C. L. and R. J. Hanson., 1974, Prentice-Hall, Inc., Englewood Cliffs., New Jersey.
Lay, T., et al., 2011, Geophys. Res. Lett, 38 (6). Art. No. L06302 . ISSN 0094-8276
Newman, A. V., et al., 2011, Geophys. Res. Lett., 38, L05302, doi: 10.1029/2010GL046498.
Satake, K., et al., 2011, to be submitted to PAGEOPH 2011.
Yagi, Y., et al., 2004, Bull. Seism. Soc. Am., 94, 1795-1807.

6