Earthquake risk assessment of the Opak and Merapi-Merbabu active faults to support mitigation program in Yogyakarta province and its vicinity

Based on historical records, Yogyakarta has a high seismic risk related to the earthquake events along active faults, such as the Opak and Merapi-Merbabu Faults. These faults were responsible for several destructive earthquakes in Yogyakarta City and its vicinity and caused fatalities and building damage in the area, e.g., the 2006 (Mw 6.3) Yogyakarta earthquake and earlier in 1943 and 1867. A previous study shows that the Opak Fault has a geodetic slip-rate of 5 mm/y and a potential magnitude Mw 6.6. In addition, the active Merapi-Merbabu Fault has a geodetic slip-rate of 1 mm/y and a potential magnitude Mw 6.6. We used scaling law relations of earthquake parameters and magnitude scenarios to estimate the recurrence time of each fault based on a kinematic model. Our results estimate that the earthquake return period (Tr) for the Opak Fault (Mw 6.6) is ∼162 years, the maximum intensity is ∼VII-VIII MMI scale, the Peak Ground Acceleration (PGA) is ∼36 % g, and Peak Ground Velocity (PGV) is ∼ 30 cm/s for a 5 km hypocentral depth. In the meantime, the earthquake return period for the Merapi-Merbabu Fault (Mw 6.6) is estimated to be ∼810 years, the maximum intensity is ∼ VI-VII MMI, the PGA is ∼ 30-36 % g, and the PGV is ∼ 21-24 cm/s for a 5 km hypocentre depth. Both faults potentially produce destructive earthquakes (Mw > 6.0) in Yogyakarta City and its vicinity. Therefore, assessments of (paleo) earthquakes are needed of both the Opak Fault and the Merapi-Merbabu Fault to support the long-term earthquake hazard mitigation program.


Introduction
Yogyakarta Province in Indonesia had been shaken by many large earthquakes sourced from the Opak Fault which is located in the eastern area of Yogyakarta City. Within several hundred years, several earthquakes destroyed infrastructure in the Yogyakarta and area around it and caused many fatalities. The destructive earthquake in 2006 (Mw 6.3, hypocentral depth 12.5 km) [1] caused more than 5,600 fatalities and around 200,000 different types of buildings collapsed. The earthquake source was located along the Opak Fault which stretches from the Southern Yogyakarta to the Northern Klaten City (Central 2 Java). The distribution of infrastructural damage had an identical direction to the strike of the Opak Fault. Another destructive earthquake occurred in 1943 with Mw 7.0, the epicenter located in the offshore of Southern Yogyakarta Province (Figure 1). That earthquake caused 213 fatalities, 2,096 injured, and around 28,000 buildings collapsed. Additionally, a destructive earthquake hit Yogyakarta City and its vicinity in 1867 ( Table 1). The exact epicenter had not been determined, but based on the historical records about 372 buildings collapsed and 5 people died [2].
There are no earthquake sources data exist in the catalog data for the last 200 years [2] for the Merapi-Merbabu Fault, but geodetic data of the Merapi-Merbabu Fault suggest a geodetic slip rate of 1 mm/y [3]. It means the Merapi-Merbabu Fault is one of the active faults in Central Java that corresponds to a potential magnitude of the earthquake of about Mw 6.6 [3]. Based on the historical seismicity and damage, we pointed out the need of a paleo seismic assessment and monitoring of Opak and Merapi-Merbabu Faults.  The tectonics of Java (Yogyakarta Region, part of Southern Java) are dominated by the subduction of the Australian Plate to the north-northeastward beneath the Sunda Plate with a relative velocity of about 5-7 cm/year [6]. The Australian Plate dips to the north-northeastward from the Java Trench, attaining depths of 100-200 km beneath the Java Island and depths of 650-700 km to the north of the island. The earthquake of 26 May 2006 at 22:53:58 UTC occurred at a shallow depth and its epicenter was estimated along the Opak Fault. The Opak Fault is located in the eastern part of Yogyakarta City. It is a strike-slip fault. It has an almost south-north direction from the Southern Coast of Yogyakarta to Prambanan Temple in the north. The neighboring fault is the Merapi-Merbabu Fault (Central Java, it is not showed in fig.1). The fault is located between Merapi and Merbabu mountains and has geodetic slip rate of about 1 mm/y [3].

Data
We used the Indonesian earthquakes catalogue data [1,2,4,5] for all of the earthquakes that occurred along the Opak Fault ( Table 1). The geodetic slip-rate data of the Opak Fault which has been published from [3] reveals about 5 mm/y and the potential earthquake magnitude is about Mw 6.6. We used this data (Mw 6.6) and converted to seismic moment to determine the static stress drop and recurrence time estimation of the earthquake in the future

Fault Length and Rupture Area Estimation
The Yogyakarta earthquakes that occurred in 1867 (Mw 8-estimated) and 2006 (Mw 6.3) indicated that the Opak Fault is active and has the potential to drive large earthquakes in the future. Therefore, the foreshocks and aftershock distribution of the 2006 earthquake (figure 1) were important data to estimate the fault length (L) and rupture area (A) in this paper. Meanwhile, to predict the Tr with targeted magnitude Mw 6.6 and less than Mw 6.6 [3], this paper used the scaling law relation of the earthquake to calculate the L and A. We assumed that the focal mechanism of the earthquake was the strike-slip model, then we estimated the L and A for Mw 6.6 or less using the Wells and Coppersmith [7] formula such as shown in equations (1) and (2).
(1) and (2) Annotation: a = -3.55; b = 0.74; c = -3.42 and d = 0.9 for strike slip fault By assuming the average value of L and A, it is possible to calculate the width of fault (W) by divided equations (2) by (1).

Static Stress Drop Theory
The tectonic energy will be released during a co-seismic event after the critical stage of rock is passed over. This energy is accumulated along the inter-seismic period since the last large earthquake event has occurred. Figure 2 shows the illustration of built-up stress and released stress. It was begun from the initial stress (σo), the build-up stress along the inter-seismic period until it achieved the critical stress (σ1) stage. After the critical stress stage passed over, it would release several stresses called a stress drop.  Figure 2. Illustration of the stress drops and return period interval of earthquake relations. The time interval between two events of earthquakes and the difference between the critical stress (σ1) value and initial stress (σ0) value is called return period (Tr) and stress drop, respectively. Modified from Kanamori [8].
Furthermore, to describe the stress release effects, we used one type only, i.e. a static stress drop. Then, we used the Kanamori and Anderson [9] formula as follows; (3) µ is rigidity (SI), is strain change, and C is a non-dimensional shape factor. To determine the shape factor (C) and L, we referred to shallow infinite transverse shear (dip-slip) faults. Kanamori and Anderson [9] referred to Starr (1928) and Aki (1966) for C = formula, where λ is the Lame constant. We used the definition of a seismic moment as; = Mo is a static seismic moment (Nm) and S is a rupture area (m 2 ). Based on the substitute equation (3) and (4), and assumed µ = λ, we found the static stress drop formula for dip slip fault as: Equation (5) is similar to Beeler and Hickman [10] which reveals that static stress drop is calculated using seismic moment measured from seismogram and estimated rupture area, while the recurrence interval is estimated from historic records.

Seismic Energy and Moment Magnitude Relation
The relation between the surface-wave magnitude (Ms) and the total energy of seismic waves (Es) had been studied by scientists such as Gutenberg and Richter (1956), Bath (1958), and several revisions Gutenberg and Richter (1956) which is revealed by Kanamori and Anderson [9]. They obtained the relation as: Log Es = 1.5 Ms + 11.8 Then, we used the conversion Mw to Ms with the formula as: Mw = 0.9239 Ms + 0.5671 (7) The equation (7) was used for magnitude interval 6.2 ≤ Ms ≤8.7 with R 2 = 0.8183.

Results and Discussion
As explained above, we used the potential maximum magnitude scenarios of Mw 6.6 for Opak and Merapi-Merbabu Faults by Pusgen (2017) [3] to calculate some fundamental of earthquake risk if the  Table 2 shows the relations between earthquake parameters, such as Mw scenario, Mo, L, and A. These parameters were used to estimate the static stress drop (Δσ) value for focal scenarios in the strike-slip model. The slip model scenario was needed to select a formula to calculate L and A based on Wells and Copper Smith [7].

Recurrence Time Estimation
We used slip rate geodetic data to estimate the earthquake recurrence time with the potential of magnitude scenarios of Mw 6.6. Davidson and Scholz [12] and Kanamori and Anderson [9] calculated the recurrence time of the earthquake by dividing the moment of an earthquake by source time function or moment rate accumulation. By using slip rate geodetic data of the Opak Fault (5 mm/y) and Merapi-Merbabu Fault (1 mm/y), the Tr could be estimated. Table 3 presents the estimation result of the return period of the earthquake that occurred in active Opak and Merapi-Merbabu Faults for several magnitude scenarios. Then diagram (Figure 4) describes the recurrence time of the active Opak Fault with slip rate based on slip rate geodetic data of 5 mm/yr [3] and other parameters based on Table 2.   Figure 5 explains the recurrence time estimation for Merapi-Merbabu active fault with using a slip rate geodetic data of 1 mm/y [3] and the other parameter.

Potential of Shaking Impact of Earthquake
One of the secondary effects of the earthquake is a shaking impact on infrastructures and the surrounding area. Figure 6 shows our models simulate the shake map, Peak Ground Acceleration (PGA), and Peak Ground Velocity (PGV) caused by strike-slip earthquake scenario (Mw 6.6) of the Opak Fault with the focus depth (h) of about 15 km. The shaking impact, PGA, and PGV were revealed in the MMI scale, in percentage of g (in % g), and in cm/s respectively. This model assumed that the hypocenter was a point source located along the Opak Fault, and it did not consider a focal mechanism. If the earthquake scenario is Mw 6.6 and depth 15 km, it will have the shaking impact of about VII MMI, PGA about 28 % g, and PGV of about 21 cm/s, respectively.      The Yogyakarta earthquake on 26 May 2006 with Mw 6.3, depth 12,5 km, and maixmum intensity of VIII MMI by USGS [1], the constructed shakemaps in figure 6, 7, 8, and 9 offered an appropriated estimation. We compared our modeled shaking impact with the field survey results conducted after the 2006 Yogyakarta (Mb 5.9 or Mw 6.3) by Murakami and Pramitasari [13]. They reported the infrastructure damage from south to northeast follows the strike orientation of the Opak Fault. The heavily damage area are spread out starting from Pundong, Imogiri, Bantul, Plered, Piyungan until the Prambanan Temple area. The damage distribution is also consistent with UNOSAT estimation published on their web-site soon after the earthquake [13]. Moreover, the damage distribution matched with the microzonation survey result conducted by Karnawati et al. [14]. They revealed that the very high susceptible zone for future earthquakes is situated in the western part of Pleret, central part of Jetis, Imogiri, and Pundong. All sites were located in Yogyakarta province. Based on this information, we suggest that our modeling result of earthquake impact shaking of the Mw 6.6, depth 5 km, and 15 km, with the epicenter was assumed as a point source at the Opak Fault agreed with the field damage impact of the 2006 Yogyakarta earthquake.
The modeled PGA with Mw 6.6 scenario yield a maximum ~36 % g(depth 5 km) and 28 % g (depth 15 km). We compared this results with field survey results from the Mw6.3 2006 Yogyakarta earthquake (Mb 5.9) conducted by Fathani (2007) [15]. He used the BMKG data (Mb 5.9) and USGS data (Mw 6.3) for 2006 the Yogyakarta earthquake to estimate PGA from several empirical formulas by Donovan (1973), Esteva (1974), Matuschka (1980), Campbell (1981), and Fukushima-Tanaka (1990). In one of his summaries from BMKG data, he determined PGA using the Esteva formula and he found the g (min) value was about 0.233 g (in BPKP site) and g (max) value was 0.322 g (Watu site). Another references were Matuschka (1980) and Fukushima-Tanaka. Matuschka (1980) found that the g(min) value was about 0.241 g and g (max) value was about 0.311 g, while Fukushima-Tanaka (1990) found that the g(min) value was 0.194 g and g(max) value was 0.287 g. These values corresponded with the PGA modeling as explained above. If compared with the earthquake in 1867, which caused the damage in Yogyakarta City, we underwent a difficulty for more interpretation because of incomplete data. The water castle at the Sultan's Palace in Yogyakarta city was heavily damaged, and Gedung Agung as Dutch Governor's residence collapsed. Visser (1922) estimated that the location of the 1967 Yogyakarta earthquake epicenter was following the Opak River [5]. Based on the damaged data, it could be estimated that the magnitude of 1867 earthquake was probably larger than Mw 6. The relations of shaking impact on MMI scale and g value were interpreted by Housner (1970), Gere and Shah (1984), and Yeats et al (1997) in Fathani (2007) [15]. Referring to their result for Ml 6.0, the earthquake impact had a scale of VII-VIII MMI and the ground acceleration was about 0.22 g. Therefore, we interpreted 10 that all shakemap results with Mw 6.6 ( Figure 6-9), which show the intensity, was about VII-VIII MMI and the g minimum was about 24 % g.
The high value of estimated MMI, PGV and PGA value suggest that the future earthquake along this fault may cause significant damageto buildings. The possibility of earthquake-induced landslides also needs to be considered to utilize early warning system, monitoring and mitigation program in the future. The Mw6.3 2006 Yogyakarta earthquake caused landslides in Prambanan District and Sengir [16]. With a comparable event in the future scenario of earthquake along the Opak Fault with more than 6.0, the hazards associated with earthquake induced landslides should be considered in some areas, primarily in the high slope area.

Conclusion
Based on historical earthquakes, Yogyakarta Region is prone of earthquake hazard region. At least two destructive earthquakes occurred and caused serious damaged in the Yogyakarta City and some areas around it, namely the 2006 earthquake (Mw 6.3) and the 1867 earthquake (unknown magnitude) which likely originated from the Opak Fault. If the active Opak Fault had a slip-rate geodetic of about 5 mm/y and estimated maximum credible earthquake about Mw 6.6, the return period (Tr) of earthquake will be around 162 years (figure 4). Given these parameters, the estimation of the maximum intensity is around VII-VIII MMI, while the PGA is around 36 % g, and PGV is around 30 cm/s (figure 7). Meanwhile, for the Merapi-Merbabu Fault which has a geodetic slip-rate of 1 mm/y and maximum probable earthquake of Mw 6.6, the estimation of earthquake return period (Tr) is around 810 years ( Figure 5). These provide an estimation of the maximum intensity value at VI-VII MMI, PGA ~ 30-36 % g and PGV ~ 21-24 cm/s ( Figure 9). Results from this study can be utilize to support the landscape planning and the long-term earthquake hazard mitigation program in Yogyakarta region and its vicinity.