Seismic Hazard Function (SHF) study of coastal sources of Sumatra Island: SHF evaluation of Padang and Bengkulu cities

Significant earthquakes on the island of Sumatra, Indonesia, have predominantly been earthquakes with a thrust mechanism that occurred due to the subduction process and seismotectonics near coastal cities of West and South Sumatra, which could be affected by earthquakes triggered by these seismic sources. We compared the Seismic Hazard Function (SHF) of two coastal cities of Sumatra: Bengkulu and Padang. The results showed that the SHF of Bengkulu is higher than that of Padang. Estimated earthquake hazards are presented in the form of seismic hazard maps expressed as the PGA of 10% rate of exceedance probability in 50 years. In estimating the seismic potential in Sumatra, the seismic moment rate was jointly estimated from the smoothed mean seismicity rate and the pre-seismic subduction surface strain rate model. In this study, the island of Sumatra was chosen as a master model for Seismic Hazard Analysis (SHA). The motivation for choosing Sumatra for the SHA was because of the large body of complete historical earthquake data of the North Western Sunda Arc. The SHF is calculated based on a magnitude range of 6.0 to 9.0 during 50 years with the radius distance from the source less than or equal to 100 km.


Introduction
Understanding seismic hazard is basically determined by how well or how reliable a potential earthquake model functions; while the reliability of a potential seismic model is determined by the extent of our knowledge regarding understanding possible sources, data completeness, and the rate of seismicity. Seismic potency in a certain area could be estimated on the basis of the seismic moment rate. The seismic moment rate can be estimated on the basis of seismic wave amplitude, GPS, or geodetic data, and the slip-rate of a Late Quaternary fault. Since the possible size of an earthquake is related to the fault length, understanding the possible segmentation is important in seismic hazard study and analysis. Seismic hazard is often estimated from seismic activity in an area of interest, using the Gutenberg-Richter magnitude-frequency relation (Gutenberg and Richter 1944). Hazard estimations obtained from this kind of approach are strongly dependent on the level of available knowledge of the seismic history of the study area and usually includes much subjective judgment.
To overcome the shortcomings of the classical approach to seismic hazard estimation, Wesnousky et al. (1983Wesnousky et al. ( , 1984 proposed the integration of geologic and paleoseismologic data. Ward (1994) proposed incorporating information from geology, paleoseismology, space geodesy, observational seismology, and synthetic seismicity. Wesnousky et al. (1984) attempted to integrate geologic and paleoseismologic data to obtain better Open Access *Correspondence: wtriyoso@gmail.com † Wahyu Triyoso, Aris Suwondo, Tedi Yudistira and David P. Sahara contributed equally to this work 1 Global Geophysics Group, Faculty of Mining and Petroleum Engineering, Bandung Institute of Technology, Jl. Ganesha 10, Bandung 40132, West Java, Indonesia Full list of author information is available at the end of the article Triyoso et al. Geosci. Lett. (2020) 7:2 estimations of the rate of earthquake occurrences. Ward (1994) tried to quantify seismic potential using geodetic data where the location and slip rates of faults are not well-known. Frankel (1995) proposed spatial smoothing of instrumentally recorded seismicity by combining small and moderate earthquakes with uniform background seismicity to avoid the subjective judgment necessary for drawing seismic source zones where causative structures of seismicity are largely unknown. Triyoso and Shimazaki (2012) proposed using a combination of the seismicity smoothing developed by Frankel (1995) and surface strain rate deduced from GPS data to estimate the potential seismic model for seismic hazard study and analysis. All proposed models are intended to improve the estimation of the seismic hazard in an area where the period of earthquake catalog observation is relatively short, GPS stations are rare, and the knowledge of Late Quaternary faults is even less. In this study, the evaluation and mapping of the Seismic Hazard Function were done by realizing the leastsquare collocation method in which GPS pre-seismic data are used as the reference for deriving a seismic moment rate model and addressed by subduction earthquake data. Furthermore, by avoiding any strain surplus and deficit, the seismic moment rate model is normalized and then weighted, using the mean seismicity smoothing rate with correlation distances of 25, 50, and 150 km. The total probability of selected sources with a radius distance ≤ 100 km from a site along the Sumatra Subduction Zone is calculated for events with a magnitude range of 6.0 to 9.0. The Seismic Hazard Function of the integrated megathrust sources model was mapped by calculating a PGA of 10% exceedance probability for 50 years. Petersen et al. (2007) constructed a seismic hazard map for Southeast Asia, including Sumatra and Java, by combining the source model's background seismicity, subduction zone segment and crustal faults. The earthquake data used are a combination of the EVC, EHB, PDE, and ISC catalogs for a greater or equal to M5 (1964-2006). The EVC is the IASPEI Centennial catalog compiled by Engdahl and Villaseñor (2002). Further developments with additional data were included in the 2017 PUSGEN updated seismic hazard map for Indonesia by implementing the same algorithm. The result of this study is comparable to previous study results for the subduction zone source (Petersen et al. 2007, The 2017PuSGeN 2017; the pattern of estimated seismic hazard seems to be similar, in which hazard frequency is high in the coast boundary area from the northwestern to the southern part of Sumatra.

Data and methods
The seismicity data used in this study were taken from the 2017 PuSGeN earthquake catalog and moment magnitude conversion. PuSGeN is an Indonesian research consortium specializing in geohazard, and consists of experts from government institutions, academic universities, and the private sector. The catalog is a compilation of several catalogs; i.e., the United States Geological Survey (USGS), the International Seismological Centre Global Earthquake Model (ISC-GEM), Engdahl, van der Hilst, and Buland (EHB), and data from the Indonesian Meteorology, Climatology, and Geophysical Agency (BMKG). The advanced double-difference relocation technique has also been applied, using regional BMKG networks to improve the accuracy of the hypocenter locations.
The GPS data used around the Sumatran islands in this study are mainly taken from Bradley et al. (2017), Chlieh et al. (2007), Prawirodirdjo et al. (2010, and Shearer and Burgmann (2010). In the case of the Andaman and Nicobar Islands, the pre-seismic velocity crustal movement model was made on the basis of forward modeling (Okada 1985(Okada , 1992 by referring to the co-seismic data of Subarya et al. (2006) and Chlieh et al. (2007) in which the plate convergence rate of 14 mm/year (Bilham et al. 2005) is used.
Earthquakes occurring with a magnitude M w ≥ 4.6 and depth range of H ≤ 50 km around the island of Sumatra dating from 1963 to 2016 were selected. We found that the regional b-value is ~ 1.0. The earthquake catalog was then declustered in order to obtain the independent earthquake events using ZMAP software (Wiemer 2001). The surface strain estimates are based on leastsquares collocation (LSC), and the result of existing GPS data around Sumatra and its surroundings (Bilham et al. 2005;Bradley et al. 2017;Chlieh et al. 2007;Subarya et al. 2006). The seismicity smoothing was also estimated for the declustered catalog data with correlation distances of 25, 50 and 150 km in which M w ≥ 5.0 and depth range of H ≤ 50 km were used for this study. Figure 1 shows the earthquake catalog data and the b-value (a) and GPS data used for surface strain rate estimation (b).

Seismicity smoothing
In this study, Frankel's algorithm (1995) of seismicity smoothing of earthquake data is used to determine A-value with minimum subjective judgment. Seismicity smoothing needs to be done because the location where the earthquake would most likely occur in the future will not necessarily be the same place where the previous earthquake occurred. Therefore, a factor which takes into account the uncertainty of future earthquake locations was used.
First of all, gridding needs to be done on the area to be studied; the number (n i ) of earthquake events with magnitudes greater than the reference (M ref ) is then counted in each cell. The count of n i represents the maximum likelihood estimate of 10 a for earthquakes above M ref in the cell (Bender 1983). The grid of n i values is then smoothed spatially by using a Gaussian function with correlation distance c. For each cell i, the smoothed value is obtained from: in which ñ i is normalized to preserve the total number of events, ∆ij is the distance between the ith and jth cells, and c is the correlation distance. In Eq. (1), the sum is taken over cell j within a distance of 3c from cell i.

Occurrence rate function
The theoretical earthquake occurrence rate function for a particular cell, be decided from the viewpoint of magnitude completeness. Thus, applying the Gaussian function to smooth the seismicity implies accepting the 10 a by Eq.
(2). Furthermore, the following equation can also be written when substituting 10 a of Eq.
(2) in Eq. (1): is the smoothed value for cell i of the number of earthquakes above reference magnitude during the time interval T, and b is the uniform b-value.

Least-squares collocation
Least-squares collocation (LSC) is a generalized estimation method that combines adjustment, filtering and prediction (Mikhail and Ackermann 1976). This method is particularly appropriate for determining the terrestrial gravity field from arbitrary data, but it can also be applied to interpolation and transformation problems that arise in geodesy. Referring to a systematic and fairly comprehensive elementary presentation of the theory and its application to the Japanese Islands (El-Fiky 1998; Oware 1998), the method of LSC was applied to define the seismic moment rate around the Sumatran Subduction Zone. In this study, we mainly adopt the approximation by Ward (1994) and others (Molnar 1979;Savage and Simpson 1997;Field et al. 1999); as such, it is applied to the LSC based surface strain rate model to calculate  (2010). In the case of the Andaman and Nicobar Islands, the pre-seismic velocity crustal movement model resulted on the basis of forward modeling (Okada 1985(Okada , 1992 by referring to the co-seismic data of Subarya et al. (2006), Chlieh et al. (2007) in which the plate convergence rate of 14 mm/year (Bilham et al. 2005) is used Triyoso et al. Geosci. Lett. (2020) 7:2 the scalar moment rate which can be expressed by the following formula: in which, µ is the rigidity, H is the seismogenic depth, A is the unit area, and e 1 and e 2 are the principal strain rates.

Hazard calculation: probability of exceedance
The annual exceedance probability of peak horizontal ground acceleration or velocity (PGA or PGV) u at a site due to events at a particular cell k under the Poisson distribution is given by: where P k (m ≥ m (u o , D k )) is the annual exceedance probability of earthquakes in kth cell, m (u o , D k ) is the magnitude in kth source cell that would produce an PGA or PGV of u o or larger at the site, and D k is the distance between the site and the source cell. In this study, D k is calculated on the basis of the distance of source to site and the top of the starting locking depth, which is a 3-km depth (The 2017 PUSGEN 2017). The function m (u o , D k ) is the Ground Motion Prediction Equation (GMPE) relation which will be discussed in the next section. The probability distribution of PGA or PGV at the site was determined by integrating the influences of the surrounding source cells, as in: By substituting the GMPE, we could obtain which gives the annual exceedance probability of particular PGA or PGV. For the specific time duration T, the probability of exceedance is given by: The annual probability of exceeding specified ground motions is calculated by applying Eq. (7) for each grid. For the specified time duration T, the probability of exceeding specified ground motions is calculated using Eq. (8).

Ground Motion Prediction Equation (GMPE)
To construct the seismic hazard map expressed by PGA or PGV, we need an attenuation relationship (Ground Motion Prediction Equation), in terms of PGA or PGV as a function of magnitude and distance. Unfortunately, (4) M o = 2µHA max(|e 1 |, |e 2 |), there is no specific GMPE that was derived for the Indonesian region. Therefore, in this study, we used GMPE that was derived for other regions or worldwide data which had similar geological and tectonic conditions and focused on megathrust, i.e., the GMPE of Fukushima and Tanaka (1992), Youngs et al. (1997) or, Zhao et al. (1997), and Atkinson and Boore (2006). Later we called the four tested GMPE as FT92, YG92, ZH97 and AT06, respectively. To select the appropriate GMPE, the SHF on each seismic cluster (presented in Fig. 2) based on four GMPE mentioned earlier were evaluated. The SHF around Padang and Bengkulu city showed similar pattern. The results of the SHF based on four GMPE around Bengkulu city are shown in Fig. 3. This graph shows that the result of the SHF based on GMPE of ZH97 and AT06 are more appropriate and reasonable compare to the SHF based on GMPE of FT92 and YG97. In which in FT92 and YG97 the PGA is saturated in the low probability of exceedance range. It might significantly affect the PGA calculation in a longer return period, i.e., 500 and 2500 years. Therefore, we used ZH97 and AT06 for PGA calculation. We could address also this selection with the database used to derive the GMPE, in which ZH97 was developed on the basis of recorded data mainly from crustal and subduction interface and in-slab earthquakes in Japan, with supplementary data from western part of the United States and 1978 Tabas, Iran, earthquakes, while AT06 was developed on the basis of real recorded ground motion data from interface and in-slab earthquakes occurring in subduction zones of Alaska, Chile, Cascadia, Japan, Mexico, Peru and the Solomon islands.

Results and discussion
The result of the LSC surface strain rate model is then used for calculating the scalar moment rate based on Eq. (4); the probability of occurrence can be estimated using Eq. (5). Finally, the seismicity rate is estimated by the normalized seismic moment rate model and then weighted, using the mean seismicity smoothing rate of M w ≥ 5.0 with correlation distances of 25, 50, and 150 km. The result can be seen in Fig. 2. The SHF curve was obtained by plotting the probability exceedance value of earthquake events versus PGA values. The total probability of selected sources with a radius distance of about 100 km along the Sumatra Subduction Zone is calculated for events with a magnitude range of 6.0 to 9.0. In the calculation of PGA, the influential parameters are the changes in the magnitude and the distance of the source, based on the GMPE of ZH97 and AT06, and calculate the mean. Thus, using Eq. 8, the PGA of 10% exceedance probability for 50 years at the base rock can be calculated (Cornel 1968;McGuire 1976;Grandori et al. 1984;WGCEP 1995). Figure 4 shows the seismic hazard maps expressed as the PGA of 10% of exceedance probability in 50 years. On the basis of the result, in terms of ground shaking, hazard is high around the coastal boundary area from the northwestern to the southern part of the island of Sumatra. The possibility of frequent high shaking occurs around the southern part of the Sumatra Islands. The Seismic Hazard Function of the nearby coastal cities of Padang and Bengkulu was then compared. Figure 5 shows that the SHF of Bengkulu is higher than Padang.

Conclusions
The estimation and analysis of the earthquake hazard function against megathrust along the island of Sumatra was carried out based on an integrated model. The Seismic Hazard Function of the integrated megathrust Seismic hazard maps expressed as the PGA of 10% of exceedance probability in 50 years. The total exceedance probability is calculated on the basis of selected sources within a radius distance of about 100 km along the Sumatran Subduction Zone for events with a magnitude range of 6.0 to 9.0. In calculating PGA, the parameters of influence are the changes in the magnitude and distance of source which are based on the GMPE of ZH97 and AT06 and the calculated mean PGA sources model is mapped by calculating PGA of 10% exceedance probability for 50 years at the base rock. In term of ground shaking, the hazard is high around the coast boundary areas from the northwestern to the southern part of the island of Sumatra. The possibility of frequent high shaking occurs around the southern part of the Sumatra Islands. The comparison result of the SHF of nearby coastal cities shows that the SHF of Bengkulu is higher than that of Padang. The results obtained in this study could be further integrated with stress heterogeneity analysis (Sahara et al. 2018) to better understand the seismic hazard potential around the off-coast and near-coastal cities of the island of Sumatra. show that the SHF of Bengkulu is higher than that of Padang