Citation: Octantyo, A., H. Miyake (2020), Strong ground motion simulation of the 2019 Java, Indonesia,earthquake (Mw 6.9) using empirical Green’s function method, Synopsis of IISEE-GRIPS master report.

STRONG GROUND MOTION SIMULATION OF
THE 2019 JAVA, INDONESIA, EARTHQUAKE (Mw 6.9)
USING EMPIRICAL GREEN’S FUNCTION METHOD
Ardian Yudhi OCTANTYO1,2Supervisor: Hiroe MIYAKE3

ABSTRACTThis study estimated strong motion generation area parameters and simulated the strong ground motionof the 2 August 2019 (Mw 6.9) intraslab earthquake and a hypothetical plate boundary earthquake (Mw8.7) using the strong motion records of the 11 August 2019 (Mw 5.1) plate boundary earthquake withthe empirical Green’s function method. We first estimated the strong motion generation area that
reproduced the strong ground motion during the Mw 6.9 intraslab earthquake, and conducted a roughestimation of the scaling parameters to see the impact of these parameters to the performance of thesynthetic waveform reproduction. We also adjusted the stress drop ratio parameter (C) for some station
records, which is sensitive to waveform amplitude. We then modeled the strong motion generation areaof the Mw 8.7 hypothetical plate boundary earthquake using the Mw 5.1 plate boundary earthquake andthe Mw 6.9 intraslab earthquake as the element events. The estimated peak ground acceleration values
from the empirical Green’s function simulation was compared with the conventional ground motion
prediction equation. The comparison provided simulated acceleration level excited by a hypothetical
Mw 8.7 that may occur in the future. The variability of simulated ground motion in terms of methods,
source type as the element event, and rupture directivity effect was confirmed.
Keywords: Simulation, Strong Ground Motion, Peak Ground Acceleration, Ground Motion Prediction
Equation, Empirical Green’s Function Method.
1. INTRODUCTION
The Indonesia archipelago is one of the regions with the most concentrated seismic activity on the Earth.
Five major tectonic plates interact within the Indonesia archipelago, in accordance with the MORVEL
model (DeMets et al., 2010). One of the highest risk regions threaten by the seismic hazard is the Javaisland, the most populated island in Indonesia. The Java island region is overshadowed by the potencies
for the next large earthquakes occurrence from the seismic activities of the off southern Java coast
subduction zone. There is a possibility of seismic gap area lead by the indication of the interplate
coupling activities (Hanifa et al., 2014) and lack of seismicity (PUSGEN, 2017) with the accumulated
seismic moment equivalent to a Mw 8.7 earthquake. It is important to understand the strong groundmotion characteristic triggered by a large earthquake, to anticipate recurrence of similar or larger
earthquakes. But since the available historical information of the strong ground motion is limited, strong
ground motion simulation may be helpful to fulfill the requirement of earthquake engineering
application as well as the seismic hazard assessment.
This study applied the empirical Green’s function method formulated by Irikura (1986) andIrikura and Kamae (1994) to perform the strong ground motion simulation for the most recent significant
earthquake event of the 2 August 2019 Java intraslab earthquake (Mw 6.9) utilizing strong motion
recordings of the 2019 Java plate boundary earthquake (Mw 5.1 and Mw 5.2) as the element events.
Earthquake source parameter modeling is accomplished by following the recipe formulated by Miyakeet al. (2003). We also tried to simulating strong ground motion for the maximum magnitude scenario
earthquake in the same region (Mw 8.7) and performing the comparison of the simulated Peak Ground1

Agency for Meteorology, Climatology and Geophysics, Republic of Indonesia.
IISEE-GRIPS master course student.
3
Associate Professor, Earthquake Research Institute, The University of Tokyo.
21

Acceleration (PGA) value obtained from the ground motion simulation with conventional groundmotion prediction equations (GMPEs).
2. DATAWe selected the 2 August 2019 (Mw 6.9) intraslab earthquake as the target event for the strong groundmotion simulation. However, this earthquake had a deficiency of aftershock and there is also noappropriate event having a similar characteristic with this earthquake in adjacent area. We decided to
select plate boundary earthquakes occurred on 11 August 2019 with magnitudes Mw 5.1 and Mw 5.2, as
the alternative element events to simulate the target earthquake, which were recorded at many strong
motion observation stations. In present study we analyze the target event which were recorded at 47
INSMN stations located in the Sumatera and Java regions as well as the couple of earthquake events
that act as the alternative Green’s function element event.

Figure 1. Hypocenters of the 2 August 2019 target earthquake along with alternative
element earthquakes (parameter and focal mechanism by BMKG, 2019) and the 47
recorded strong motion stations (yellow dots).
3. THEORY AND METHODOLOGYThis study applied the empirical Green’s function method by the formulation of Irikura (1986) andIrikura and Kamae (1994) to obtain well-performed strong ground motion simulation. This method
utilized smaller event records to reproduce the ground motion simulation of a larger event, with theassumption that those small event records consolidate propagation paths and local site effects (Poiataand Miyake, 2017). The scaling relations between the target event and the element one in the empiricalGreen’s function formulation are described into the ratio of fault dimensions (N) and stress drop (C)parameters, which is important to determined beforehand to obtain simulation of strong ground motion.
The parameters of N and C are determined by identifying the flat-level at the displacement
and acceleration spectra amplitude of the two events. By the source spectral fitting, the estimation of themoment ratio of the target event to the element one event (M0 / m0), corner frequencies of the target
event (fcm) and the element one (fca) are provided. The relationships among the obtained parameters withthe scaling parameters of N and C are described in Equations (1) and (2) as follows.
𝑀0= 𝐶𝑁 3
𝑚0𝑀𝑓(f → 0)2(𝑚0 ) ( 𝑓𝑐𝑚 ) = 𝐶𝑁0𝑓

𝑁 = 𝑓𝑐𝑎

𝑐𝑚𝑐𝑎and(f → ∞)𝑀𝑓(1)3𝐶 = (𝑚0 ) ( 𝑓𝑐𝑚)
0𝑐𝑎(2)The target event waveform is reproduced by summing the waveforms of element eventswith adjustment for the difference in the slip velocity time function between the element event and thetarget event. The equation of summation is written as Equation (3),2𝑁
𝑈(𝑡) = ∑𝑁
𝑖=1 ∑𝑗=1

𝐹(𝑡) = 𝛿(𝑡 − 𝑡𝑗 ) +1
1𝑒𝑛′(1− )𝑟𝐹(𝑡) ∗ (𝐶 ∙ 𝑢(𝑡))𝑟𝑖𝑗1∑(𝑁−1)𝑛′[𝑘=1(𝑘−1)𝑒 (𝑁−1)𝑛′𝑟 −𝑟𝑜𝑡𝑖𝑗 = 𝑖𝑗

𝑉𝑠(𝑘−1)𝑇𝛿 {𝑡 − 𝑡𝑖𝑗 − (𝑁−1)𝑛′}]

𝜉+ 𝑖𝑗

(3)

𝑉𝑟

where U(t) represents the synthetic waveform of the target event, while u(t) the observed waveform ofthe element event, r the epicentral distance from each recording stations, Vs refers to S-wave velocity
around source area and Vr refers to rupture velocity along the fault plane, while T is the rise time of thetarget event, N and C are the ratio of the fault dimension and stress drop of the target event to the element
one, respectively, and F(t) is the filtering (adjustment) function for the difference in the slip velocity
time function.
4. RESULTS AND DISCUSSION
4.1. Source parameters estimation for the 2 August 2019 (Mw 6.9) earthquake
Waveform data of the target event, as well as the waveform data of the 11 August 2019 (Mw 5.1)
earthquake and the 11 August 2019 (Mw 5.2) earthquake as the element events, from the CLJO andKASI stations, were used to estimate the scaling parameter of fault dimension ratio (N) and stress drop
ratio (C). We conducted a rough estimation of the scaling parameters to see the impact of these parameterto the performance of the synthetic waveform reproduction. We found that the stress drop ratio (C)parameter is sensitive to the amplitude of the synthetic waveform estimation for this pair of earthquakes.
We had to assign the stress drop ratio (C) parameter separately for KASI station located in the azimuth
direction of 350o and other stations that located in the azimuth direction range of 41 o – 91o so thesimulation of the synthetic waveform can be well performed. Detailed adjusted scaling parameters value
is summarized in Table 1.
Table 1. Scaling parameters of the 2/08/2019 (Mw 6.9) target
event pairs with the 11/08/2019 (Mw 5.1) element event.
otherNo.
scaling parameters
KASIstations
1fcm
0.3 Hz0.3 Hz2fca

0.9 Hz0.9 Hz3N

334

C

10

105

propagation pathadjustment coefficient
best fit C

1

3

1030

6

Figure 2. Strong motion generation
area estimated to simulate the Mw
6.9 intraslab earthquake using theMw 5.1 plate boundary earthquake.

Figure 2 shows the schematic model of the strong motion generation area for the target
event. Sub fault dimension was set to be 2.0 km in length by 2.0 km in width, with the rise time of 0.2
second. Therefore, with the fault dimension ratio (N) equal to 3, the strong motion generation area
dimension of the target event became 6.0 km in length by 6.0 km in width, with the rise time of 0.6
second and the initial rupture starting point at sub fault at the shallowest northern end of the strong
motion generation area (1, 1).
Comparisons of the observed and synthetic waveforms for the KASI station and the CLJO
station are provided in Figure 3. Here we found that through the adjustment application of the stressdrop ratio (C), adequate synthetized waveform of the estimated strong motion generation area can beachieved for acceleration, velocity and displacement pulses recorded at the two stations.

3

The adjustment for the parameter of stress drop ratio (C) was done to prevent too large missestimation of the synthetic waveform reproduced from the element event waveform, since we used plateboundary earthquake records as the element to synthetize the waveform of the intraslab earthquake.
The difficulty to apply only one single source parameter as the input for the ground motion
simulation performance can be influenced by the factors of lateral variations, attenuation structure, orother heterogeneity of subsurface structure along the propagation path, since we have significant
different in depth, fault distance and azimuth degree between the pair of earthquakes. Another factor
that might be involved is the differences in the source mechanism that are assumed to have homogeneousradiation pattern due to the scale of the magnitude.

Figure 3. Ground motion simulation results of the CLJO and KASI stations (EW component).
4.2. Source parameters estimation for the Mw 8.7 hypothetical plate boundary earthquake
We tried several simulations scenario for the hypothetical plate boundary earthquake with Mw 8.7. First
simulation utilizing the plate boundary earthquake Mw 5.1 as the element event, with the closest sub
fault from the element event hypocenter as the rupture starting point (7,1). We also tried to assign theopposite direction rupture starting point (7,5) to see the sensitivity due to rupture propagation directivity.
In the first step of the first simulation, we simulated a temporary earthquake event of Mw 6.9; wespecified the N parameter of 10 and C parameter of 1. In the second step of the first simulation, weutilized the synthetic waveform from the temporary earthquake event, as the element input for thesimulation of the target hypothetical plate boundary earthquake of Mw 8.7. We specified the final N
parameter of 7 in length and 5 in width with C parameter of 1.
Determination of N parameter was done using scaling relations and other supporting
information from previous studies. For C parameter we were consistently using the value of 1 byassuming same stress drop value between target and element events, and consider the simulations as thelower bound simulations. The C parameter controlled the amplitude of the simulated waveform, when
we assigned larger C parameter, we will expect larger simulated peak ground acceleration.
Second simulation utilizing the intraplate earthquake Mw 6.9 as the element event, to findthe possibility of using large intraslab earthquakes, that often occur inside the deeper part of the plate,to simulate the larger plate boundary earthquakes. Since we wanted to compare the results with the firstsimulation, we assigned one of the rupture starting points that is used in the first simulation, as thestarting point of the second simulation. We also used the N and C parameters that already determined in
the first simulation, without any adjustment. Radiation pattern correction was applied for the secondsimulation. The schematic illustration for the strong motion generation area of the Mw 8.7 hypotheticalplate boundary earthquake are shown in Figures 4 and 5.
Simulated acceleration waveform and peak ground acceleration values of the hypothetical
Mw 8.7 plate boundary earthquake for CLJO and KASI stations shows in Figure 6. The simulation result
derived from the starting point (7,1) generally has larger PGA value compared to those derived from thestarting point (7,5), that indicating the rupture directivity effect to the ground motions in each station.

4

Figure 5. Strong motion generation area
estimated to simulate the Mw 8.7 hypotheticalplate boundary earthquake.

Figure 4. The Mw 8.7 hypothetical plateboundary earthquake simulation scheme.

56.244.2

36.2

CLJO

CLJO

CLJO

181.6

171.7

KASI

KASI

76.0

KASI

Figure 6. Simulated acceleration waveforms of the hypothetical Mw 8.7 plate boundary earthquake at
CLJO and KASI stations (NS component). 1st column is the results for the 1st simulation with rupture
starting point (7,1), 2nd column is for rupture starting point (7,5), and 3rd column is the results from the
2nd simulation with the rupture starting point (7,5).
4.3. The application of simulated
ground motion: long-period groundmotion in Jakarta
We applied the simulation of the Mw
8.7 hypothetical plate boundary
earthquake using the Mw 6.9 intraslab
earthquake as the element event in
Jakarta. Here we would like to find thelong-period ground motion contents
and identify the characteristic
behavior in the synthetic waveform
Figure 7. Comparison of the pseudo velocity response in
resulting from simulation process. In
Jakarta station between the Mw 6.9 intraslab earthquake (leftFigure
7, we can notice from the flat
figure) and the simulated Mw 8.7 hypothetical plate boundary
levelofthe pseudo velocity responseearthquake (right figure).
that the Mw 6.9 intraslab earthquake
already contains long-period ground motion in Jakarta station ranging 0.5 - 7.0 seconds, and it wasenhanced in the application result of the simulation for the Mw 8.7 plate boundary earthquake.
5

4.4. Comparison of the simulated PGA values from strong ground motion simulation and groundmotion prediction equation (GMPE) for the Mw 8.7 hypothetical plate boundary earthquake
Comparisons between the results
obtained from empirical Green’sfunction simulation with groundmotion prediction equation (GMPE)
by Youngs et al. (1997) for the Mw 8.7
hypothetical plate boundary earthquake in CLJO and KASI stations,showed that GMPE calculation in
CLJO station for rock condition was
more comparable to the empiricalFigure 8. Comparison of the estimated PGA values
Green’s function simulation. Whileobtained from the empirical Green’s function simulations
for KASI station, GMPE calculation
and GMPE by Youngs et al. (1997) for plate boundary
for soil condition was more
source in rock and soil site conditions.
comparable to the empirical Green’sfunction simulation. These differences were influenced by the local site condition in each station.
Detailed estimated PGA values obtained from various methods are presented in Figure 8.
5. RECOMMENDATION
Hopefully in the future, the study of GMPE for Indonesia archipelago, can be well established alongwith the accumulation of strong ground motion records. So that, we can have better understanding on
the strong motion attenuation relationship. Then, we can discuss the comparison with the empiricalGreen’s function simulation from the point of view of PGA values.
We confirmed the variability of simulated ground motion in terms of various
methodologies that we used to estimate the PGA values in this study as well as the variation of thesource types that we used as the element event, and also the rupture directivity effect.
GMPE method and the empirical Green’s function method often give different results,
which bring us to a question “which one we should use as the appropriate estimation of PGA?”. In thiscase, GMPE is more dependent, controlled by the constant, parameters and standard deviation that wasalready provided by the equation. On the other hand, empirical Green’s function simulation is more
independent; we can assign arbitrary parameters to obtain the simulation of the ground motion, andsometimes it becomes dangerous. So that, N and C parameters are crucial as the simulation control
parameters. Supporting information regarding the parameters from previous study prior to the simulation
will be helpful. Other idea is conducting trial simulations for all rupture starting points and calculate theaverage PGA value along with the standard deviation, to better fit the seismic hazard assessment.
ACKNOWLEDGEMENTS
This research was conducted during the individual study period of the training course “Seismology,Earthquake Engineering and Tsunami Disaster Mitigation” by the Building Research Institute, JICA,and GRIPS. I would like to acknowledge my sincere gratitude to my supervisor Dr. Hiroe Miyake andmy advisor Dr. Toshiaki Yokoi for their perpetual guidance and support during this research work.
REFERENCES
DeMets, C., Gordon, R. G., and Argus D.F., 2010, Geophys. J. Int., 181, 1-80.
Hanifa, NR., Sagiya, T., Kimata, F., Efendi, J., Abidin, HZ., and Meilano, I., 2014, Earth Planet. Sci. Lett.,
401, 159-171.
Irikura, K., 1986, Proceedings of the 7th Japan Earthquake Engineering Symposium, 151-156.
Irikura, K., and Kamae, K., 1994, Ann. Geofis., 37, 1721-1743.
Miyake, H., Iwata, T., and Irikura, K., 2003, Bull. Seis. Soc. Am., 93, 2531-2545.
National Center for Earthquake Studies of Indonesia (PUSGEN), 2017, Ministry of Public Works andHousing, Republic of Indonesia.
Poiata, N., and Miyake, H., 2017, Pure Appl. Geophys., 174, 3503-3519.
Youngs, R. R., Chiou S.-J., Silva W. J., and Humphrey J. R., 1997, Seismol. Res. Lett., 68, 58-73.

6