Citation: Simborio, T., T. Yokoi, T. Hara (2022), Strong ground motion simulation of the 2020 Masbate,
Philippines earthquake (Mw6.6) using Empirical Green’s Function method, Synopsis of IISEE-GRIPS
Master’s Thesis.

STRONG GROUND MOTION SIMULATION OF THE
2020 MASBATE, PHILIPPINES EARTHQUAKE (Mw6.6)
USING EMPIRICAL GREEN’S FUNCTION METHODTom Carlo E. Simborio1
MEE21705

Supervisor: Toshiaki YOKOI2 , Tatsuhiko HARA3
ABSTRACT

Strong ground motions of the August 18, 2020 Mw 6.6 Masbate, Philippines earthquake were simulatedusing the August 11, 2020 Mw 4.9 Masbate, Philippines earthquake records of the Philippine Strong
Motion Network through the Empirical Green’s Function (EGF) method. The fault dimensions andstress drop ratios were calculated through the source spectral fitting method. A grid search was done toobtain a set of parameters for the EGF method that explains the observed strong motion data well. These
parameters include the rupture starting point, rupture velocity, rise time, length, and width of the strong
motion generation area (SMGA). The result suggests the rupture started in the shallow southern segmentof the fault plane. The peak ground accelerations (PGAs) from the simulated waveforms are consistent
with the observed PGAs, and a directivity effect was also simulated. The PGAs calculated using aground motion prediction equation used for hazard estimations in the Philippines overestimates for threestations among the four stations analyzed in this study. The peak ground displacements (PGDs) fromthe simulated waveforms are underestimates compared to the observed PGDs. The ratio of the size of
the SMGA to the moment magnitude and the ratio of the rise time to the moment magnitude are smaller
than those obtained in previous studies.
Keywords: Strong Ground Motion, Empirical Green’s Function Method, Peak Ground Acceleration.
1. INTRODUCTION
The Philippines experiences earthquakes every day due to its geotectonic setting. Within the country is
the Philippine Fault Zone (PFZ) that transects the entire country and divided into segments. One of its
segments is the Masbate segment. Significant earthquakes occurred in the history of Masbate; onedamaging earthquake occurred on August 18, 2020, with a moment magnitude Mw 6.6.
As of November 2021, the Philippines has 103 strong-motion stations throughout thecountry. Recorded events are analyzed, cataloged, and stored in a database for scientific research,especially in earthquake engineering, to update the National Structural Code in the Philippines (Grutas
et al., 2017).
The objective of this study is to investigate the applicability of the Empirical Green’s
Function (EGF) method to the current strong ground motion data of the Philippine Strong Motion
Network (PSMNet). Specifically, I conducted the following: I estimated the parameters for the strong
motion generation area and compared them to those obtained by previous studies. I compared the peakground accelerations (PGAs) from the EGF simulations to the observed PGAs and compared the peakground displacements (PGDs) from the EGF simulations to the observed PGDs. I also compared the
PGAs from the EGF simulations to the computed PGAs based on the ground motion prediction equation(GMPE) by Fukushima and Tanaka (1990, 1992).
2. METHODOLOGY
1

Department of Science and Technology-Philippine Institute of Volcanology and Seismology (DOSTPHIVOLCS), Philippines
2
Japan International Cooperation Agency (JICA)
3International Institute of Seismology and Earthquake Engineering, Building Research Institute (IISEE-BRI)

1

2.1. Empirical Green’s Function (EGF) method
Empirical Green’s Function (EGF) is a method to synthesize waveforms of a large event, which isreferred to as a target event, using seismograms recorded from a smaller event, which is referred to as
an element event (e.g., Irikura, 1986; Irikura et al.; 1997; Miyake et al., 2003). In the EGF method, thefault plane of the target event is divided into NxN subfaults, and the subfault corresponds to the size ofthe fault plane of the element event. The waveform of a target event, U(t), is calculated by:N
U(t) = ∑N
i=1 ∑j=1

r
F(t) ∗ [C ∙ u(t)]
rij

(1)

where U(t) is the synthetic waveform of the target event; i and j are indices of the subfaults; r is the
hypocentral distance from the element event to the recording stations; rij is the distance from recording
stations to each subfaults; F(t) is the Modified Filtering (Correction) function; C is the ratio of the stressdrops between the element and target events, and u(t) is the observed waveform of the element event.
The Modified Filtering function F(t) by Irikura et al. (1997) is used.
2.2. Source spectral fitting method
As stated above, the parameters N and C are necessary to perform the EGF method. Miyake et al. (1999)
proposed the source spectral fitting method to determine N and C parameters objectively. This methodfits the average of the observed spectral ratio between the target and element events to a theoreticalMsource spectral ratio to estimate ( o ), corner frequencies of the target event (fcm) and element eventmo(fa), where Mo and mo are the seismic moments of the target and element events, respectively.
For calculations of observed source spectra, I set Vs to 3.2 km/s and used Qs(f) =105.02f0.75 which Oestar and Hayashida (2021) obtained. I minimized the sum of the squared
differences between the average of the observed spectral ratio and to a theoretical source spectral ratio
over a frequency range, while Miyake et al. (1999) took the standard deviations of the observed spectral
ratio into account. The minimization is conducted by a heuristic method introduced by Press et al.
(1992), which combined the Down Hill Simplex Method (Nelder and Mead, 1965) and the Very Fast
Simulated Annealing (Ingber, 1989). N and C are calculated by:f

N = f ca

(2)

cm

M

3

f

C = (mo ) ( fcm )
o

(3)

ca

2.3. Grid search
After the values of N and C are obtained, the grid search is carried out to find parameters of the EGF
to explain the observed strong ground motion records well. This method is done by initially setting
parameters for the position of rupture starting point, rupture velocity, a rise time, the length, and widthof the subfaults; executing EGF, and calculating residuals. The misfit function to be minimized is:
Residual = ∑Station ∑Component

2
2
∑N
∑N
i=1(Uobs − Usyn )i=1(Aenv. obs − Aenv. syn )
+5·[
2
2]
2
N
2
N
2
2
(∑N
(∑N
i=1 Uobs · ∑i=1 Usyn )i=1 Aenv. obs · ∑i=1 Aenv. syn )

6

(4)

{
}
where Uobs and Usyn are the observed and synthetic displacement, respectively; Aenv. obs. and A env. syn. are
the envelopes of the observed and synthetic acceleration, respectively.
I put a larger weight for the differences of the acceleration envelopes, because thedifferences of the displacements are larger. I determined the weight in Eq. (4) after trying several
different weights.
3. DATA
3.1. Earthquake parameters2

Two earthquakes are necessary to perform the Empirical Green’s Function (EGF) method: a target event(mainshock) and an element event – a smaller event that occurred in the vicinity of the mainshock.
Through inspection, I selected the August 11, 2020 event (“08/11/2020 – 0726Z” in Figure 1), which isnearest to the mainshock (“08/18/2020 – 0003Z” in Figure 1) among the earthquake events as theelement event. The distance between the foreshock and mainshock is approximately 4 km. The
earthquake source parameters of the target and element events are summarized in Table 1. The second
nodal plane (NP2) was used in this study based on the surface rupture and aftershock distribution.
Table 1. The origin times, hypocenters, magnitudes, and the focal mechanisms of the element
and target events. The origin times and hypocenters are from the earthquake catalog of the DOSTPHIVOLCS. The magnitudes and focal mechanisms are from the moment tensors solutions fromthe SWIFT system (https://swift1.phivolcs.dost.gov.ph/).
Date

Time (UTC) Latitude (N°) Longitude (E°) Depth (km) Magnitude (Mw) NP

Element
08/11/2020 07:26:41.23
Event

11.943124.055

13.2

4.9

TargetEvent

11.969124.027

11.4

6.6

08/18/2020 00:03:46.79

Figure 1. Moment tensor solutions of the2020 Cataingan, Masbate earthquake,and its foreshocks and aftershocks.

Strike (°) Dip (°) Rake (°)

NP1

231

85

177

NP2
NP1

321
235

87
64

5
-174

NP2

142

85

-26

Figure 2. Station distribution ofaccelerometers that recorded theforeshock and mainshock.

3.2. Strong motion records of the foreshock and mainshock and data preparation
Twenty Kinemetrics ETNA-2 accelerometers of the PSMNet recorded both the foreshock andmainshock. I used four stations since data from at least four stations are necessary to perform EGF
calculations (Miyake et al., 2003). I chose the closest stations in Masbate, Bicol Region, Leyte Island,
and Cebu Province (MMPA, BCPL, LYKN, and CLIL, respectively) from the events (Figure 2). The
recorded waveforms are saved as an Event Log file. These records are converted to ASCII files using
the ViewWave software (Kashima, 2022) and are converted into files in the KNET format(https://www.kyoshin.bosai.go.jp/kyoshin/man/knetform_en.html). I added appropriate headers to the
KNET format files and obtained ZOO format files. Data in the ZOO format are necessary to perform
Empirical Green’s Function method.
I used a time series starting from 8 seconds before the theoretical P-arrival time, whose
duration is 81.92 seconds. I manually picked an S-arrival time using a strong motion record and its Husid
plot (Husid, 1969). A cosine taper function in which the value of the function becomes one at the Sarrival time is multiplied to the strong motion data. The tapered part is extracted and used in this study.
I used the vector sum, which is the square root of the sum of the squares of the acceleration
amplitude spectra of the horizontal components. Parzen window was used to smoothen the acceleration
amplitude spectra to reduce the oscillatory behavior of the spectra.

3

4. RESULTS AND DISCUSSION
4.1. Source spectral fitting
Figure 3 shows the result of thesource spectral ratio fitting. The
thick theoretical curve shows thefrequency range used for fitting. I
used the average of the source
spectra ratios in the frequency
range from 0.2 to 8 Hz for source
spectra ratio fitting. The lowfrequency asymptote is 0.01,
whereas the high-frequency
asymptote is 0.2. The cornerfrequency of the mainshock is
0.33 Hz, whereas the cornerfrequency of the foreshock is 1.33
Hz. From Eq. (2) and (3), it can becalculated that the ratio of fault
dimensions (N) is 4, and the ratioof the stress drop between theevents (C) is 1.56.

Figure 3. The result of
the source spectral ratio
fitting.

Figure 4. Schematic illustrationof the strong motion generation
area of the August 18, 2020,
Masbate, Philippines earthquake.
The star represents the rupture
starting point.

4.2. Grid search
I performed a grid search for parameters of the EGF method. The parameters searched in this study are
the position of rupture starting point, rupture velocity, a rise time, and the length and width of thesubfaults. Based on the model parameters that minimize the misfit function defined in Eq. (4), the strikeof the rupture is well-constrained at position 3, the dip is at position 1, the rupture velocity is 3.0 km/s,
the rise time is 0.05 s, length, and width of the subfaults are 1.1 km. Figure 4 shows a schematic
illustration of the strong motion generation area of the August 18, 2020, Masbate, Philippines earthquake
where the star denotes the rupture starting point. The result suggests that the rupture started at the shallow
southern segment of the fault plane.
4.3 Empirical Green’s Function simulation, peak ground accelerations (PGAs), and peak grounddisplacements (PGDs)
Synthetic seismograms of the target event were computed by the EGF method using the set of parameters
obtained by the grid search method. An example is shown in Figure 5. The observed and synthetic
acceleration waveforms are similar to each other. The differences between the displacement waveforms
are large.
Figure 6 shows the comparison of the different PGAs. I compared the observed PGAs
with the synthetic PGAs and PGAs obtained from the ground motion prediction equation (GMPE) by
Fukushima and Tanaka (1990, 1992). I also compared the observed PGDs with the synthetic PGDs.
Comparing the observed PGAs to the synthetic PGAs, they are comparable. Comparing the observed
PGA to the PGA computed by the GMPE, the latter overestimate for three of the four stations analyzed
in this study. In terms of PGDs (Figure 7), the synthesized PGDs tend to underestimate.

4

Figure 5. The observed and synthetic waveforms of LYKN (N-S)

Figure 6. Comparison of the observed PGAs and Figure 7. Comparison of the observed PGDs andthe calculated PGAs.
calculated PGDs.
4.4 Source scaling
I compared the results of this study to the empirical source scaling relationships for crustal earthquakes
by Sommerville et al. (1999) and Miyake et al. (2003). The comparison is made by comparing the
moment magnitude with the SMGA size (Figure 8 (a)) and with the rise time (Figure 8 (b)). For the
moment magnitude of 6.6 of the 2020 Masbate earthquake, the SMGA size and rise time are smaller
than the general trends of the previous studies.

Figure 8. Comparison of the moment magnitude to (a) SMGA and (b) rise time.
5. CONCLUSIONS

5

This study has shown the applicability of the Empirical Green’s Function (EGF) method to the current
strong ground motion data of the Philippine Strong Motion Network (PSMNet). The parameters for thestrong motion generation area (SMGA) were estimated through grid search. The result suggests that therupture started in the southern segment of the fault plane, and the rupture is shallow. The peak groundaccelerations (PGAs) from the EGF simulations to the observed PGAs were comparable, and theobserved effect of the directivity is reproduced. The PGAs from the EGF simulations were compared tothe computed PGAs based on the ground motion prediction equation (GMPE) by Fukushima and Tanaka(1990, 1992), and it shows the PGAs calculated using the aforementioned GMPE tend to overestimate.
The peak ground displacements (PGDs) from the EGF simulations to the observed PGDs tend tounderestimate. Compared to previous studies, it shows the size of the SMGA and rise time is smaller.
Further applications will help examine the source characteristics of earthquakes in and around thePhilippines.
ACKNOWLEDGEMENTS
This research was conducted as an individual study for the training course “Seismology, Earthquake
Engineering, and Tsunami Disaster Mitigation” by IISEE-BRI, JICA, and GRIPS. I sincerely thank Dr.
Toshiaki Yokoi and Dr. Tatsuhiko Hara for supervising me. I thank Dr. Irikura and Dr. Miyake for theirlectures on Strong Ground Motion Study and the Empirical Green’s Function program. I also thank the
DOST-PHIVOLCS for nominating me for this opportunity and assisting in my data requests. I used the
Seismic Analysis Code (SAC, Goldstein et al., 2003; Goldstein and Snoke, 2005) in some parts of thedata processing in this study. The Generic Mapping Tools (GMT, Wessel and Smith, 1998) was used tomake a figure in this study.
REFERENCES
Goldstein, P., Dodge, D., Firpo, M., Minner, L., 2003, edited by Lee, WHK, Kanamori, H., Jennings,
P.C., and Kisslinger, C., Academic Press, London.
Goldstein, P., Snoke, A., 2005, Incorporated Institutions for Seismology Data Management Center
Electronic Newsletter.
Grutas, R., Tiglao, R., Lasala, M., Deocampo, J., Narag, I., and Solidum Jr., R.U., 2017, IAG-IASPEI.
Fukushima, Y., and Tanaka, T., 1990, Bull. Seism. Soc. Am., 80, 757-783.
Fukushima, Y., and Tanaka, T., 1992, The Seismological Society of Japan 1992 Fall Meeting.
Husid, L.R., 1969, Características de Terremotos, Análisis General, Revista del IDIEM 8, Santiago,
Chile, 21–42.
Ingber, L., 1989, Mathematical and Computer Modelling, 12, 967-973
Irikura, K., 1986, Proceedings of the 7th Japan earthquake engineering symposium, 151–156.
Irikura, K., Kagawa, T., and Sekiguchi, H., 1997, Seism. Soc. Japan 2, B25 (in Japanese).
Kashima, T., 2022, ViewWave (version 2.2.6), http://smo.kenken.go.jp/~kashima/viewwave.
Miyake, H., Iwata, T., & Irikura, K., 1999, Zisin, 51, 431–442. (in Japanese with English abstract).
Miyake, H., Iwata, T., and Irikura, K., 2003, Bull. Seis. Soc. Am., 93, 2531 – 2545.
Nelder, J. A., and Mead R., 1965, The Computer Journal, 7, 4, 308–313
Oestar, J., and T. Hayashida, 2021, Synopsis of IISEE-GRIPS Master’s Thesis.
Press, W.H., Teukolsky, S.A., Vetterling, W.T., and Flannery, B.P., 1992, Cambridge University Press,
USA.
Somerville, P., Irikura, K., Graves, R., Sawada, S., Wald, D., Abrahamson, N., Iwasaki Y., Kagawa T.,
Smith, N., Kowada, A., 1999, Seismological Research Letters, 70, 59–80.
Wessel, P., and Smith, W.H.F., 1998, New, improved version of Generic Mapping Tools released, Eos,Trans. AGU, 79, 579, 1998

6