Citation:
Rosid, A., Y. Yagi (2021), Complex behavior in the source process of the 18 august 2020,
Bengkulu, southwest of Sumatra doublet earthquake, Synopsis of IISEE-GRIPS Masterʼs Thesis.

COMPLEX BEHAVIOR IN THE SOURCE PROCESSOF THE 18 AUGUST 2020,
BENGKULU, SOUTHWEST OF SUMATRA
DOUBLET EARTHQUAKE
Abdul Rosid1,2

Supervisor: Yuji YAGI3
ABSTRACTThis study investigated the source process of the 2020 Bengkulu doublet earthquake (Mw 6.9 and Mw
7.2 by this study) using the high degree of freedom finite fault inversion method recently developed by
Shimizu et al. (2020). In the process, we applied some unique treatments by considering an identicalfault plane and using the same data set for both earthquakes. Consequently, we have succeeded inreconstructing a reliable interaction of the foreshock and mainshock source processes by showing good
waveform fitting between the synthetic waves and the observed waves.
We found that the rupture in the fault plane was divided into two distinct areas between
foreshock and mainshock. The foreshock earthquake broke the asperity located in the southeast with a
total maximum slip of 1.07 meters, while in the mainshock earthquake, the asperity broke with a total
maximum slip of 1.17 meters in the northwest direction. Furthermore, the extension pattern of the 2020
Bengkulu doublet earthquake rupture shows as if the mainshock earthquake continued the unfinishedforeshock earthquake rupture but with a more comprehensive area coverage. It seems that the final stageof the foreshock rupture triggered the mainshock hypocenter, which is closely located. Meanwhile, thesuspending pattern of the foreshock's and mainshock's rupture extension to the deeper part appears to berelated to the 2007 Bengkulu Great Earthquake (Mw 8.5). Moreover, using the high degree of freedom
finite fault inversion method, which allows us to extract the fault plane geometry information, we
discover that the dip angle changes along the fault plane as the depth increases.
Keywords: doublet earthquake, source process, 2020 Bengkulu earthquake, rupture interaction, trigger.
1. INTRODUCTION
Indonesia is an archipelago country that has quite complicated tectonic conditions. This condition resultsfrom a convergence of the world's major plates and several small plates or micro blocks (Bird, 2003).
As a result of these conditions, Indonesia also has a high seismic activity every year. Recently, anearthquake with a moment magnitude (Mw) 6.7 (foreshock) occurred in the early morning on 18 August2020 in the Southwest of Sumatra. According to the Global CMT (https://www.globalcmt.org), theorigin time of this earthquake is 22:24:12 UTC and located in 4.56 S, 100.99 E, with a depth of 23 km.
Soon after that, another earthquake occurred in almost the exact location with almost the same
magnitude. It was Mw 7.0 (mainshock) earthquake that occurred around 5 minutes after the firstearthquake. It was located at 4.47 S, 100.86 E, with a depth of 25 km and its origin time is 22:29:38
UTC. By considering the parameters of the two earthquakes, it can be said that this is a twin or doublet
earthquake phenomenon if we refer to Lay and Kanamori (1980). With the doublet earthquake on 18
1

Meteorological, Climatological, and Geophysical Agency (BMKG), Indonesia.
IISEE-GRIPS Master’s course student.
3
Professor, Department of Life and Environmental Sciences, University of Tsukuba.
2

1August 2020, southwest of Sumatra, we will use far-distance seismic broadband records from stations
worldwide to conduct a finite fault inversion representing fault deformation on an assumed fault by
shear slip vectors by superposition of five basis double couple components to those earthquakes. Thismethod was a new method introduced by Shimizu et al. (2020). In addition, by using the same method,
we can extract the geometry information of the subduction slab along the fault plane, which is the sourceof the doublet earthquake.
2. DATA
In conducting this research, we used teleseismic body-wave data recorded by stations located aroundthe world. Those were downloaded from the official website of Incorporated Research Institutions for
Seismology (IRIS; https://ds.iris.edu/wilber3/find_event). By applying several constraints to theteleseismic waveform data, such as the distance range of the stations is 30o - 90o from the epicenter(yellow star) (Figure 1), they have a good signal-to-noise ratio in the vertical component, and they musthave good azimuthal coverage. We got 54 teleseismic waveforms data recorded from 54 stations
worldwide (red triangles) (Figure 1), both for the foreshock and the mainshock.
We deliberately utilized the same dataset for both
earthquakes, because considering the fact that both are twin
earthquakes which, of course, occurred at almost the exact
location and time. Therefore, by downloading the waveforms
starting from 1 minute before the P-wave phase of the foreshock
to 10 minutes after, we have obtained not only the seismic wave
recording data of the foreshock but also the mainshock in only one
dataset. The mainshock occurred just 5 minutes after theforeshock. In the condition of a doublet earthquake recording, theP-wave arrival of the second earthquake (mainshock) is hard to
identify since it is contaminated by the foreshock wave recording,which is still not over. Consequently, we applied a time-shifted
method in determining the P-wave arrival of the mainshock in thisstudy. Once we finish picking the foreshock, we calculate a time
Figure 1. Azimuthal equidistancedifference between the foreshock and mainshock origin time.
map view of teleseismic stations
Then, we employed this time difference to assign the P-wave
which is recording Foreshock andarrival of the mainshock automatically. This method is not only
Mainshock.
helpful in solving picking accuracy problems but is also a crucial
step in understanding the interaction of a doublet earthquake since we want to infer the mainshock
earthquake source parameters based on the foreshock earthquake picking. Furthermore, we referred tothe hypocenter and the fault parameter of the doublet earthquake from the Global Centroid Moment
Tensor (GCMT; https://www.globalcmt.org), and we assumed the velocity structure based on the 1Dmodel proposed by Collings et al. (2012) since it has the exact location of targeted earthquakes.
3. METHODOLOGY
In investigating the source processes of the doublet earthquake, we implemented a finite-fault inversion
that is more flexible than conventional inversion designs. Conventional inverse solutions have beensustained by restricting the model space and reducing the degree of freedom for slip vectors.
Nevertheless, these restrictions are not fundamentally physical elements for representing source
processes. Furthermore, improper conjectures about the fault geometry can increase modeling errors,
produce nonunique final solutions, and make the solutions obtained becomes arduous to interpret(Shimizu et al., 2020). In general, the source of an earthquake can be expressed as the volume momentrate tensor density Ṁ(t,ξ) (Backus and Mulcahy, 1976), implied by a linear combination of five basis
double couple components Mq (Kikuchi and Kanamori, 1991) as shown by Eq. (1):(1)

2Ṁ 𝑡, 𝝃ṁ𝑡, 𝝃 𝑴𝒒

where ṁq(t,ξ) is a spatiotemporal moment-rate volume-density function of qth basis component of
moment tensor. While, qth here corresponds to the degree of freedom of the moment rate tensor without
isotropic expansion sources, and ξ symbolized a location in the source area.
In considering the inversion framework of Shimizu et al. (2020), we built a synthetic
teleseismic waveform in the inversion process based on Green's function and the assumed source processmodel. Furthermore, we compared the synthetic waveform with the observed teleseismic waveform data
for each worldwide recording station to justify that the source process model obtained by this study is
relevant to its original state. A small error could indicate this. Eq. (2) shows how we built synthetic
teleseismic waveforms for individual stations worldwide regarding the observed waveform data.
𝑢 𝑡

𝑎

𝐺

𝑡, 𝜉

𝛿𝐺

𝑡, 𝜉

∗ 𝑋 𝜉 𝑇 𝑡

𝑡𝑘 𝑑𝜉

𝑒

𝑡

(2)

where uj represents seismic waveform for far-field term observed at a station j by a linear combinationof moment-rate volume-density functions of five basis components qth of double couple moment tensor,
K and L denote the slip direction on the model plane S, aqkl are the expansion coefficients to be estimatedfrom the observed data, 𝐺 𝑡, 𝜉𝛿𝐺 𝑡, 𝜉 is Green’s function, 𝑋 and 𝑇 are the basis
functions for space and time, respectively, and 𝑒 is a background and instrumental noise.
Furthermore, to perform a stable result, we applied a weighted moment tensor finite fault
inversion method proposed by Yamashita et al. (2021). The smoothness constraints formulation used in
Shimizu et al. (2020) is adjusting an identical Gaussian distribution for all basis components base on
Akaike’s Bayesian information criterion (Akaike, 1980). It possibly can bias the solution since thepotency-rate density functions of basis components have a chance to experience a disproportionately
smoothing compared to their amplitude. Potency density itself is referred to as a product of the averageslip on a fault and the fault area, and potency-rate density is time variation of the potency density.
However, by using the new method developed by Yamashita et al. (2021), we can adjust the standard
deviation of the smoothing constraints proportionally with its magnitude for each basis double-couple
component.
4. RESULTS AND DISCUSSION
We assessed the spatiotemporal distribution of the potency density tensor for the 2020 Bengkulu doublet
earthquake by implementing the flexible finite-fault inversion method to teleseismic P waveforms. We
employed a unique treatment in the inversion process by setting the same fault parameter with a strikeof 321o and dip of 11o for the foreshock and mainshock by considering them as a doublet earthquake.
Figure 2c and Figure 3a show that the fault has a length of 91 km along the strike and a width of 49 km
along the dip direction. However, the rupture of the foreshock mainly occurred only in a limited area
(red square) with a length of 49 km along the strike and 35 km along the dip direction. Furthermore, we
found that there is only one asperity broke for each earthquake. The asperity of the foreshock is a little
bit deeper and narrower than the mainshock, with the total maximum slip reached is 1.07 m. Meanwhile,
the asperity of the mainshock is more comprehensive, with the total maximum slip reached is 1.17 m.
The total moment tensor for both foreshock (Figure 2a) and mainshock (Figure 3c) shows a thrust fault
mechanism corresponding to the subduction zone activity in southwestern Sumatra.

3

a)

d)

c)

b)

Figure 2. Inversion Result of the Foreshock. a) total moment tensor. b) source
time function. c) total slip distribution. d) spatial-temporal rupture snapshot.

a)b)

Figure 3. Inversion Result of the Mainshock. a) total slip distribution.
b) spatial-temporal rupture snapshot. c) total moment tensor. d) source
time function.

c)d)Moreover, by using the high degree of freedom finite fault
inversion method, which allows us to extract the fault plane geometry
information, we found that the dip angle along the fault plane changes
with depth. If we look more closely at the shape of the focal mechanism from the shallow to the deeppart of the fault plane, we can find that the dip angle increases with increasing depth. At the shallowest
edge with a depth of 7.3 km, the initial dip angle starts from 7 degrees; then at a depth of 8.6 km, thedip angle becomes 9 degrees; at a depth of 10 km, the dip angle is 13 degrees; at a depth of 11.3 km,
the dip angle is 20 degrees, respectively. As for the deeper parts, we cannot consider the dip angle
because the slip generated in that area is relatively small.
In this study, we also set an exact hypocenter location for both foreshock and mainshock
because we want to infer the new location of the mainshock hypocenter by referring to the location of
the foreshock hypocenter considering that those are doublet earthquakes. The foreshock starts the
rupture sequence of the doublet earthquake. The rupture duration of the foreshock is 15 seconds (Figure2b), and it is started from the hypocenter (Figure 2d; black star). As the maximum rupture speed is set
at 3.7 km/s, the foreshock rupture starts to expand circularly around the hypocenter and extends to thedowndip direction. However, the propagation towards the deeper part is suspended at 5-6 seconds. Then,the rupture only extends in the shallower part in the southeastward direction from the hypocenter andreaches the peak moment release at 9 seconds (Figure 2b). Finally, the rupture ends at 15 seconds near
the hypocenter (Figure 2d).
The mainshock then continues the doublet earthquake rupture sequence. In this study, we
set the maximum rupture velocity of the mainshock as infinity because we want to infer the new locationof the mainshock hypocenter. By doing this, all areas within the fault plane can break from the beginning.
Then, we can identify the new location of the mainshock hypocenter by considering the high slip rate
which occurs in the early seconds. Figure 3b shows that the high slip rate occurs at 1-2 seconds in the4west-northwestward from the foreshock hypocenter (black star). Therefore, we determined that such alocation is the new mainshock hypocenter (blue star). If we pay attention, the mainshock hypocenter isestablished close to the termination area of the foreshock rupture. Consequently, we suppose that theforeshock rupture triggers the mainshock rupture. The mainshock rupture propagates northwestward
until 3 seconds; then, it extends swiftly toward the deeper part until it is suspended at 6 seconds at theexact area where foreshock rupture was also suspended. After that, it only expands in the shallower partand reaches the peak moment release at 9 seconds (Figure 3d). Eventually, starting from the 13th second(Figure 3b), the series of rupture processes remains limited only in the southeastern part with consistentshrinkage until the 20th second, which is the expiration course of the mainshock earthquake source
process (Figure 3d).

Thi
Figure 4. Map view of the doublet
earthquake.

Figure 5. Slip distribution model of the 2007 Bengkulu
earthquake (Mw 8.5) derived from the coseismic offset
modelling at SuGAr GPS sites (modified from
Figure 4 shows the location of the Ambikapathy et al., 2010).
foreshock hypocenter (orange star) andmainshock hypocenter (green star), including
the 2020 Bengkulu doublet earthquake rupture sequence where orange polygons denote foreshock
rupture snapshots and blue polygons represent mainshock rupture snapshot. If we look at the doublet
earthquake fault plane scene by considering the bathymetry height, we can say that the doublet
earthquake rupture occurred in the Sumatra forearc zone, specifically at the boundary between theaccretion wedge (bluish contour) and the forearc island (brownish contour). Collings et al. (2012)
discovered that the area beneath the forearc island has a different structure to the accretionary wedge
closer to the trench. It is indicated by a higher Vp/Vs ratio under the forearc island, implying that the
sediments in this area are more compressed and contain much less water which is squeezed during
compaction. (Hyndman and Peacock, 2003).
Furthermore, Figure 5 is a slip distribution model for the 2007 Bengkulu earthquake (Mw8.5) inferred from GPS data (Ambikapathy et al., 2010). we can find that the rupture of the doublet
earthquake (blue rectangle) is located just southwestward of the 2007 Bengkulu earthquake asperity(reddish rectangle) location. We suspect that the structure discontinuity between the accretionary wedge
and the forearc island is playing a role and becomes the reason why the foreshock and mainshock
earthquake ruptures suspend when trying to extend toward deeper parts. Indeed, the structure contrastof the two regions also corresponds to the difference in rock stress capacity; hence that the 2020
Bengkulu doublet earthquake ruptures that occurred in a lower stress area was inadequate to trigger abreak in the asperity area of the 2007 Bengkulu earthquake (Mw 8.5) which had higher stress capacity
where the energy was already released during the 2007 Bengkulu earthquake.

55. CONCLUSIONSThis study investigated the source process of the 2020 Bengkulu doublet earthquake using the highdegree of freedom finite fault inversion method recently developed by Shimizu et al. (2020) to theteleseismic bodywave data. We also applied several unique treatments in the inversion process,considering that the doublet earthquake has an adjacent location, time, and magnitude.
Accordingly, we found that the rupture in the fault plane was divided into two distinct areas
between foreshock and mainshock indicated by the opposing location of two asperities. Furthermore,
the movement pattern of the 2020 Bengkulu doublet earthquake rupture shows as if the mainshock
earthquake continued the unfinished foreshock earthquake rupture but with a more comprehensive area
coverage. In addition, we also suspect that the existence of a structure discontinuity located in the
deepest part of the fault plane affects the cessation of the doublet earthquake rupture movement in thatarea. Moreover, using the high degree of freedom finite fault inversion method, which allows us to
extract the fault plane geometry information, we discover that the dip angle changes along the fault plane
as the depth increases.
ACKNOWLEDGEMENTS
This research was conducted as the individual study of the training course “Seismology, Earthquake
Engineering and Tsunami Disaster Mitigation” by IISEE/BRI, JICA, and GRIPS. I want to express mysincere gratitude to my supervisor Prof. Yuji Yagi for his valuable support and experience, which helped
me to complete my research successfully. I would also like to acknowledge Dr. Bunichiro Shibazaki forhis permanent support from the beginning to the end of the program.
REFERENCES
Akaike, H. (1980). Likelihood and the Bayes procedure. Trabajos de Estadistica Y de Investigacion
Operativa. https://doi.org/10.1007/BF02888350
Ambikapathy, A., Catherine, J. K., Gahalaut, V. K., Narsaiah, M., Bansal, A., and Mahesh, P. (2010).
The 2007 Bengkulu earthquake, its rupture model and implications for seismic hazard. Journal of
Earth System Science. https://doi.org/10.1007/s12040-010-0037-2
Backus, G., & Mulcahy, M. (1976). Moment tensors and other phenomenological descriptions ofseismic sources—II. Discontinuous displacements. Geophysical Journal of the Royal Astronomical
Society. https://doi.org/10.1111/j.1365-246X.1976.tb01275.x
Bird, P. (2003). An updated digital model of plate boundaries. Geochemistry, Geophysics, Geosystems.
https://doi.org/10.1029/2001GC000252
Kikuchi, M., & Kanamori, H. (1991). Inversion of complex body waves - III. Bulletin - SeismologicalSociety of America.
Lay, T., and Kanamori, H. (1980). Earthquake doublets in the Solomon Islands. Physics of the Earthand Planetary Interiors. https://doi.org/10.1016/0031-9201(80)90134-X
Okuwaki, R., Hirano, S., Yagi, Y., and Shimizu, K. (2020). Inchworm-like source evolution through ageometrically complex fault fueled persistent supershear rupture during the 2018 Palu Indonesia
earthquake. Earth and Planetary Science Letters. https://doi.org/10.1016/j.epsl.2020.116449
Shimizu, K., Yagi, Y., Okuwaki, R., and Fukahata, Y. (2020). Development of an inversion method to
extract information on fault geometry from teleseismic data. Geophysical Journal International.
https://doi.org/10.1093/gji/ggz496
Yamanaka, Y., and Kikuchi, M. (2003). Source process of the recurrent Tokachi-oki earthquake on
September 26, 2003 inferred from teleseismic body waves. Earth, Planets and Space, 55(12).
https://doi.org/10.1186/BF03352479
Yamashita, S., Yagi, Y., Okuwaki, R., Shimizu, K., Agata, R., and Fukahata, Y. (2021). Consecutive
ruptures on a complex conjugate fault system during the 2018 Gulf of Alaska earthquake. Scientific
Reports. https://doi.org/10.1038/s41598-021-85522-w

6