Seismic Hazard Function Mapping Using Estimated Horizontal Crustal Strain Off West Coast Northern Sumatra

A seismic hazard study and analysis of the megathrust source off the west coast of North Sumatra, Indonesia, were conducted based on the estimated horizontal crustal strain using the surface displacement data. This area was selected due to the availability of pre- and co-seismic Global Positioning System (GPS) data for the 2005 Nias–Simeulue Mw 8.6 event. This study aimed to estimate the seismic hazard function (SHF), which is expressed as peak ground acceleration (PGA) versus probability of exceedance (PE), for a 500 years return period using GPS data. The source area model of the Mw 8.6 event is determined based on the co-seismic GPS data. The horizontal crustal strain of the source area is estimated using least square prediction employing local covariance functions based on the horizontal displacement data. The Mw 8.6 return period is estimated by dividing the sum of the co-seismic seismic moment by the pre-seismic seismic moment based on GPS data. The seismicity rate model above a magnitude of completeness is then estimated assuming the b-value of 1 obtained on the previous study’s earthquake catalog data in the region. We show that the SHF based on the study area’s horizontal crustal strain is higher than the one based on earthquake catalogs and estimated geological sliprate data. This discrepancy is associated with the static stress increase (Coulomb failure stress, CFS) of about 0.25 bar imparted by the 2004 Aceh Mw 9.1 event that occurred in the north of the study region. We interpreted that the increase of the SHF was due to the increase in the region’s stress load, which was well documented by the GPS data.


INTRODUCTION
It has been concluded that the Sumatra Subduction Zone is one of the most active plate tectonic margins in the world (Petersen et al., 2004). Its geometry has an oblique NE-ward convergence between the subducting Indian-Australian Plate and the overriding southeastern Eurasian Plate (McCaffrey, 1991;Prawirodirdjo et al., 2000;Pan et al., 2001). A convergence rate of about 49 mm/year (Zachariasen et al., 2000) results in Sumatra Island having a very high annual rate of earthquakes: over the last 250 years, five major earthquakes (M w ≥ 8.0) have occurred along the Sumatran megathrust (Megawati and Pan, 2009). The major earthquakes might perturb the vicinity stress concentration and, hence, the seismic hazard function (Steacy et al., 2005;Pollitz et al., 2006).
Previous regional hazard models for slip along the Sumatran megathrust have been proposed. The Global Seismic Hazard Assessment Program (GSHAP., 1999) proposes that two sources, i.e., the Sumatra Subduction Zone and the Sumatran Fault, characterize this region. Megawati and Pan (2001) further expand this ground motion assessment toward Singapore and Malaysia. Petersen et al. (2004Petersen et al. ( , 2007 updated the GSHAP map (GSHAP., 1999) by compiling the updated earthquake catalogs to develop new seismotectonic models. Bilham et al. (2005) use geological data to infer the seismic hazard in this region. However, none of these studies integrate geologic, seismic, and geodetic data to estimate the seismic hazard function of Sumatra Island. Ward (1994) introduced integrating information based on geology, paleoseismology, space geodesy, and observational seismology. According to Ward (1994), geology and paleoseismology come into play as they locate, segmentize, and fix the long-term slip of principal faults. Input from geodesy in the form of strain rate data takes command in quantifying the seismic potency of an area where there is inadequate knowledge of the slip rates of Late Quaternary fault data (Ward, 1998). Synthetic seismicity contributes by estimating the statistical possibilities of various scenarios of multiple segments breaking. In essence, Ward's work tries to consider an unidentified zone by using geodetic data. The essential part of his work is to incorporate geodetic data into a probabilistic seismic hazard study and analysis using Kostrov's formula regarding seismic moment and strain (Kostrov, 1974). The use of Global Positioning System (GPS) data for seismic hazard study and analysis has been widely explained in previous studies (e.g., Ward, 1994; Working Group California Earthquake Probabilities [WGCEP], 1995;Savage and Simpson, 1997). Triyoso et al. (2020) have estimated the seismic hazard function in southern Sumatra based on integrated pre-seismic GPS data, earthquake catalog, and estimated geological sliprate data. The horizontal crustal strain was estimated using the least squares collocation (LSC) technique employing local covariance functions based on the surface displacement data. The seismic moment rate model was derived from the GPS displacement of pre-seismic data around the subduction zone. To avoid the possibility of any strain surplus or deficit, the seismic moment rate model was then normalized. The rate was then used to weight the mean seismicity smoothing rate derived based on the correlation distances of 25, 50, and 150 km (Frankel, 1995), and the seismic hazard function (SHF) was constructed.
The 2005 Nias-Simeulue M w 8.6 event was recorded by the available pre-and co-seismic GPS data. Using the recorded co-seismic GPS data and Kostrov's formula (Kostrov, 1974), the source area model of the M w 8.6 event can be estimated. The return period is calculated by dividing the sum of seismic moments of the co-seismic GPS data by the preseismic GPS data for the same source area model. Furthermore, using b-value ∼1 (Triyoso et al., 2020), the annual A-value (Gutenberg and Richter, 1944) above magnitude completeness could be estimated.
Thus, this study intended to perform the probabilistic seismic hazard analysis using the source area's surface displacement data off the west coast of northern Sumatra. The seismic source model considers pre-seismic and post-seismic moment to constrain an M w 8.6 event's recurrence period. The study area is around the subduction zone of M w 8.6. Based on the previous results (Triyoso et al., 2021), the Sumatran Fault's impact could be neglected.
The seismic hazard function obtained in this study is then compared with the seismic hazard functions obtained from each single data, i.e., geology, geodetic GPS data, and seismology. The result shows that the SHF based on the horizontal crustal strain is higher than the SHF obtained from data based on earthquake catalogs and geological data. The increase of probability of exceedance (PE) in the SHF, which could be correlated with the time advance of the M w 8.6 occurrence, is further associated with the increased Coulomb failure stress in the source area of the M w 8.6 event by about 0.25 bar caused by the 2004 Aceh M w 9.1 event. This study's critical finding is that increasing static stress could increase the probability of exceedance of the peak ground acceleration (PGA) level in SHF compared to the previous study's expected result. The result may be very beneficial for earthquake mitigation and modeling efforts for probabilistic seismic hazard study and future analysis.

Data
In this study, two primary data were used to calculate the seismic hazard function of the west coast of northern Sumatra: surface displacement data based on the GPS of both pre-seismic and co-seismic and the earthquake catalog. Figure 1A summarizes the GPS and earthquake data used in this study. Following a previous study (Triyoso et al., 2020), the horizontal crustal strain values estimated based on the surface displacement data of Chlieh et al. (2007); Prawirodirdjo et al. (2010), Shearer and Burgmann (2010), and Bradley et al. (2017) were used along with additional data obtained from Bird (2003) and Qiu et al. (2019). The velocity vectors for the pre-seismic data obtained from Bird (2003) are based on GPS surveys from 1991 through 2001. Thus, the GPS data could be grouped into pre-and co-seismic data of the M w 8.6 event. Twenty particulars of the pre-seismic and 22 particulars of the co-seismic GPS data distributed in northern Sumatra were used. The horizontal displacement in each cell was then estimated using the least squares collocation (LSC) technique (Kaula, 1963;Moritz, 1978;Kato et al., 1998;Triyoso and Shimazaki, 2012;Triyoso et al., 2020). Figure 1B shows that the earthquake catalog data used in this study are based on Triyoso et al. (2020), in which seismic data and moment magnitude conversion are taken from the 2017 PuSGeN (Tim PusatStudiGempa Nasional-2017 PuSGen], 2017). The period of the seismic data is from 1963 to 2016. Following Triyoso et al. (2020), the regional b-value in northern Sumatra is assumed to be 1.0.

Least Squares Prediction Technique
Following Triyoso et al. (2020), the scalar moment rate based on Ward (1994) could be rewritten as: The µ is the rigidity, H is the seismogenic thickness, A is the unit area of the study, and e 1 and e 2 are the principal strain rates. Following Molnar (1979) and Field et al. (1999), the Gutenberg-Richter (GR) relationship in terms of the seismic moment could be expressed as: where β is equal to (b/c), b is the b-value of the GR relationship (Gutenberg and Richter, 1944), and c is the constant of the seismic moment-magnitude relationship, where log M o = cM+ d. M o could be replaced by M w using the Hanks and Kanamori (1979)relationship.
The study area is gridded into a cell size in which the surface strain rate in each cell needs to be estimated based on GPS data. We adopt procedures from previous studies to estimate each cell's horizontal crustal strain rate (Kato et al., 1998;Triyoso and Shimazaki, 2012;Triyoso et al., 2020), applying the LSC method.
As a follow-up to the successful result of the previous study in the generalized estimation of interpolation, which combines adjustment, filtering, and prediction, i.e., Kaula (1963); Mikhail and Ackermann (1976), and Moritz (1978Moritz ( , 1980, the LSC method is briefly discussed in this study. The detailed application for surface strain data is based on a previous study done around the Japanese Islands (Kato et al., 1998;Triyoso and Shimazaki, 2012) and the Sumatra Islands (Triyoso et al., 2020). The basic assumption of the LSC application is that data vector L could be composed of systematic errors (El-fiky and , in which the tectonic signal and noise could be expressed as: Following Moritz (1980), the parameter model that underlies the least square collocation is: Variable x is the vector observation or measurement, X is the unknown vector, and n is noise or the vector of the measuring errors; A is a known rectangular matrix. Thus, x could be decomposed into a systematic part denoted by AX and a random part represented by n. In the sequel, the number of observations will consistently be denoted by q and the number of parameters by m; to get an over determined problem, we must have m<q.
x and n are the column vectors of q components, X is a column vector of m components, and A is a q by m matrix. By admitting n is the noise, a second random quantity s, which is then called the signal, the following generalization could be expressed as: in which the measurement x consists of a systematic part, AX, and two random parts, s and n. Usually, the systematic part will be nonlinear initially, and the linear form AX will be obtained on linearization using Taylor's theorem. The signal s may exist at points between the measuring points, and it may vary continuously, although x is only measured at discrete points. Thus, by applying this algorithm, we could do the socalled interpolation.

Horizontal Crustal Strain Estimation
In this study, we assume that the horizontal displacement field at the observation point over the entire seismogenic thickness is homogeneous and isotropic, i.e., considering the horizontal displacement components u or v as the signal t in the equation, in which components u and v are the displacement components in the x and y directions. The local x-y coordinate system is the x-axis in the E-W direction and the y-axis in the N-S direction. A further assumption is then made in which signals u and v are not correlated with one another (El-fiky and . Even though they are correlated, the cross-covariance is very small compared to the variances (El-fiky and ; therefore, the above non-correlative assumption would be acceptable. Thus, by following the non-correlative assumption, it is possible to treat the displacement components u and v separately, thus simplifying the procedures in both the estimation of covariance function and the derivation of the least prediction equation squares collocation. By knowing the covariance functions for both displacement components u and v, we may have (El-fiky, 1998): U and V are signals to be estimated. C Uu is the (n × 1) vector of the cross-covariance between U and each element u i of u and C V v is the n*1 vector of cross-covariance between V and each element v i of v. C uu is the variance matrix of u with the size of (n × n) C vv is the variance matrix of v with the size of (n × n) u and v are the (n × 1) vectors of the known displacement components in the x and y directions. Using the notation of We may write Eq. (6) as: When u and v are determined, a and b can be obtained immediately. They are constants to U and V and, thus, can be computed beforehand. The algorithms for sparse and symmetrical linear equations are still feasible. By adopting the assumption that the covariance function is satisfied by the Gaussian function (El-fiky and  with d as the correlation distance, we may have: Then, we have: where C Uui and C V vi are the elements of the vector. The distance (d i ) between (x,y) and (x i ,y i ), with (x,y) being the coordinates of signal U or V, could be estimated by using d i = [(x-x i ) 2 +(yy i ) 2 ] 1/2 , and (x i ,y i ) are the coordinates of observation point i.
Referring to Eq. (10), we can easily obtain the following equation just by differentiating Eq. (8): where a i is the element of vector a(a 1 ,a 2 ,. . .a i . . .a n ) T and b i is the Following El-fiky and , the dilatation , maximum shear strain e max , principal strain (e 1 , e 2 ), the orientation of the major principal axis of strain θ (angle to x-axis counter clockwise), and rotation ω can be deduced from Eq. (11). Using Eqs. (11) and (1), we can calculate the seismic moment in each of our cell maps over the entire study area. We could then calculate the seismic moment of the co-and pre-seismic data based on the surface strain data.

Seismicity Rate Modeling: Source Area and Earthquake Rate Formulation
The earthquake occurrence rate above or equal to magnitude completeness as the magnitude reference (M c ) in the particular cell i could theoretically be expressed as: N i is the number of earthquakes with magnitude ≥ M c in cell i, and T is the length of record or estimated return period; v i represents the likelihood of the 10 a of the earthquake with a magnitude greater than or equal to M c (Bender, 1983). Furthermore, substituting 10 a of Eq. (12) in frequencymagnitude of the Guttenberg-Richter equation (Gutenberg and Richter, 1944), we may write the following equation: in which n i (≥M ref ) is the estimated number of earthquakes above or equal to the completeness moment magnitude. T is a given period of observation and b is the uniform b-value.

Seismic Hazard Function Estimation: Ground Motion Prediction Equation and Probability Exceedance
The PE of the annual earthquake rate of given magnitude completeness, which could be converted into the estimated ground motion [PGA or peak ground velocity (PGV)] using ground motion prediction equation (GMPE), denoted by a at a site or point of observation due to events at a particular cell k under the Poisson distribution, could be expressed as: in which P k (m ≥ m(a o , R k )) is the annual PE of earthquakes in the kth cell, m(a o , R k ) is the magnitude in the kth source cell that would produce an estimated PGA or PGV of a o or larger at the site, and R k is the distance between the site and the source cell. Since the purpose of this study was to make a comparison with previous results, the calculation of SHF, the parameter is based on Triyoso et al. (2020), in which R k is the distance of the source to the site with the starting locking depth as the top being a 3 km depth, as taken from the 2017; Tim PusatStudiGempa Nasional-2017 (The 2017 PuSGen) (2017). The function m(a o , R k ) is the GMPE relation. The total PE distribution of PGA or PGV at the site was determined by integrating the influences of the surrounding source cells, which could be expressed as: By substituting the GMPE in Eq. (15), we could obtain the annual PE of the particular PGA or PGV as follows: For a given specified time duration T, the PE could be estimated as follows: Thus, the annual PE of the specified ground motion for each grid is calculated using Eq. (16). For a given time duration T, the PE of the specified ground motions is computed using Eq. (17).

RESULTS AND DISCUSSION
This study's motivation was to perform a probabilistic seismic hazard analysis using the source area's surface displacement data off the west coast of northern Sumatra on March 28, 2005. Based on the previous study by Triyoso et al. (2021), the Sumatran Fault's impact on the seismic hazard analysis in the vicinity of the major subduction event could be neglected. Considering the availability of the recorded pre-and co-seismic GPS data in the FIGURE 3 | Estimated source are a based on the horizontal crustal strain estimated of the 2005 Nias-Simeulue event (M w 8.6) co-seismic slip (A) and the selected pre-seismic moment rate the same as that of the source area event (B). The blue dash line is the boundary area. Using the M w 8.6 event's obtained source area, we could estimate the pre-seismic moment rate's sum-up prior to M w 8.6. Knowing the seismic moment in the pre-and co-seismic data of the 2005 Nias-Simeulue, the return period (T) of the M w 8.6 event could then be estimated.  (Frankel, 1995) weighted by normalizing the estimated seismic moment rate model (C). FIGURE 5 | Detailed flow of the calculation of the mean peak ground acceleration (PGA) of 10% of exceedance probability in 50 years based on the ground motion prediction equations (GMPEs) obtained by Zhao et al. (1997) and Atkinson and Boore(2006) using the input of the seismic rate model of M w ≥ 5.0 ( Figure 4B). It shows that the epicenter of M w 8.6 took place at around the highest area of the PGA estimated. study area, we intend to compare the seismic hazard function estimated using a horizontal crustal strain with the one obtained previouslyusing various input data. A lesson learned from the comparison might be beneficial for seismic mitigation purposes.
In the previous study, the SHF of probabilistic seismic hazard analysis was constructed based on pre-seismic GPS data, data from the earthquake catalog, and the estimated geological sliprate in northern Sumatra. The geological sliprate of Bilham et al. (2005) is used to develop the pre-seismic GPS model on the position of the co-seismic GPS data of Subarya et al. (2006) and Chlieh et al. (2007). In doing so, the Okada dislocation model (Okada, 1985(Okada, , 1992 is implemented. The LSC technique using local covariance functions (Eq. 9) based on the surface displacement data is applied to estimate the horizontal displacement in each cell. The obtained GPS pre-seismic model is used to calculate the horizontal crustal strain rate and the seismic moment rate model. Afterward, the seismicity rate model is included in the algorithm to avoid the possibility of any strain surplus and deficit. It was done by normalizing the estimated seismic moment rate model and weighting the mean seismicity smoothing rate derived based on the correlation distances of 25, 50, and 150 km (Frankel, 1995). Seismicity smoothing is applied for the declustered catalog data with magnitude completeness of M w ≥ 5.0 and a depth maximum of 50 km. The declustering is performed using ZMAP (Wiemer, 2001). By applying the seismicity smoothing algorithm, we assume that the rate model is equivalent to the earthquake rate model.
As a consequence of the seismicity smoothing (Frankel, 1995) applied in this study, SHF was constructed based on point source approximation. It is applied to simplify the calculation as the realization of hazard calculation depends only on the source's magnitude and distance (Cornel, 1968;McGuire, 1976;Grandori et al., 1984;Working Group California Earthquake Probabilities [WGCEP], 1995). The source position is placed in the middle of the seismogenic thickness.
In this study, to realize the SHF of probabilistic seismic hazard analysis, we follow the GR law (Gutenberg and Richter, 1944). In this case, the annual A-value, i.e., the earthquake occurrence rate function (v i ) above the magnitude reference (M ref ), has to be estimated. M ref can be taken to be equal to or larger than magnitude completeness (M c ). Following Triyoso et al. (2020), the M ref is set at 5.0. Assuming the b-value is equal to 1, the A-value of M ref 5.0 can be estimated using the larger magnitude event rate data. Thus, the co-and pre-seismic GPS data of the FIGURE 6 | The seismic hazard function (SHF) based on the horizontal crustal strain data model (SHF-GPS) compared to the SHF based on the earthquake catalog and geological data. The SHF-GPS tends to be higher than the SHF derived from the earthquake catalog and geological data. In comparing the probabilistic seismic hazard analysis with the previous study, we pick the area with a relatively high PGA estimated in the present and the previous study and plot the SHF.
2005 Nias-Simeulue event (M w 8.6) can be used to estimate the seismicity rate by estimating the return period.
The LSC in this study is applied to estimate the horizontal displacement in every10 × 10 km 2 cell of the study area. The cell's position in this study area is the same as in the previous study (Triyoso et al., 2020). Before estimating the horizontal displacement in each cell, we need to evaluate the horizontal displacement prediction in each original data position. Figure 2 shows the results of the horizontal displacements of the co-seismic and pre-seismic GPS data. It appears that the results of the horizontal displacement estimation of both coand pre-seismic GPS are good enough, in which the predicted horizontal displacement is close to the real data. An important finding based on the LSC shows that the area of relatively high horizontal crustal strain rate took place at around the relatively high horizontal crustal strain release caused by the earthquake event on March 28, 2005. It happened around Nias Island. This phenomenon explains that the tsunami is relatively small during this M w 8.6 earthquake, or almost no tsunami was reported.
In estimating the return period, the 2005 Nias-Simeulue event's source are a needs to be defined first by summing up the seismic moment of the co-seismic GPS data. The result is then converted to the seismic moment's entire source area into M w , equal to 8.6 (Hanks and Kanamori, 1979). The model of the source area of the 2005 Nias-Simeulue is presented in Figure 3A.
Using the source area of M w 8.6, we can sum up the seismic moment of the pre-seismic GPS data, shown in Figure 3B. The return period is then calculated by dividing the sum of seismic moments of the co-seismic GPS data by the sum of seismic moments of the pre-seismic GPS data for the same source area model. The sums of seismic moments of the co-seismic and preseismic GPS data over the entire source boundary are about 118.2 × 10 27 and 0.94 × 10 27 dyn.cm, respectively. We obtained the return period of M w 8.6 as about 126 years.
Thus, using the b-value ∼1, as shown in Figure 4A, we could estimate the 10 a of M w ≥ 5.0 of the study areas as illustrated in Figure 4B. Figure 4C shows the 10 a of M w ≥ 5.0 in the previous study (Triyoso et al., 2020). Figure 4C is derived based on integrating the seismicity rate model of the mean seismicity smoothing rate derived based on the correlation distances of 25, 50, and 150 km (Frankel, 1995) weighted by normalizing the estimated seismic moment rate model.
We could then estimate the horizontal crustal strain using Eq. (11). Furthermore, the seismic moment in each cell can be estimated using Eq. (1). We assumed the rigidity (µ) and the seismogenic depth (H) to be 3.4 × 10 11 dyn.cm −2 and 20 km, respectively (Triyoso et al., 2020). The M w 's error obtained from the M w 8.6 event, based on the horizontal crustal strain estimation, could be evaluated by calculating the estimated displacement and the average difference of the observed data. The range values are small, i.e., about ∼0.1 for the co-seismic data and ∼0.2 for the pre-seismic data.
Based on the result of 10 a with M w ≥ 5.0 and the b-value ∼1, the seismic hazard map expressed as the PGA of a 500 years return period at the base rock could be constructed. The SHF curve of the total probability of an exceedance value of 10% of earthquake events in 50 years versus the PGA values is also estimated. In this study's probabilistic seismic hazard calculation, the maximum radius distance of about 100 km is used with the magnitude range of 6.0-8.6. The same magnitude range is also used to compare the SHF based on the previous study's seismicity rate model. Thus, by using Eq. (17), the calculation of PGA is realized. Following the previous study, the SHF is constructed based on point source approximation.
The GMPE used in this study refers to the recommendation results of Triyoso et al. (2020). They are the GMPEs of Zhao et al. (1997) and Atkinson and Boore (2006). It was done as we wanted to analyze the impact of the different data sources in building the SHF off the west coast of northern Sumatra. Thus, we implement both models and calculate the mean of both solutions. Figure 5 shows an illustration of the processing flow in calculating the mean PGA of 10% PE in 50 years based on the GMPE (Triyoso et al., 2020) with the input of seismic rate modeling based on Figure 4B. Figure 5 shows that the epicenter of M w 8.6 took place around the highest area of the PGA estimated. The relatively high PGA estimated in the present study is almost similar to that of the previous study, except for the closure of the highest area of the PGA estimated. The previous study shows that the PGA's highest area closure took place around the southeast of Nias Island.
In comparing the probabilistic seismic hazard analysis with the previous study, we pick the area with a relatively high PGA estimated of the present and the previous study and plot the SHF as shown in Figure 6. The previous study's earthquake seismicity rate model is based on Triyoso et al. (2020), and the geological rate is based on Bilham et al. (2005). The result shows that the SHF based on the horizontal crustal strain model tends to be higher than the SHF derived from the earthquake catalog and geological rate.
The previous study of the Coulomb failure stress (CFS) (Pollitz et al., 2006) estimated that the M w 9.1 event caused increased static stress loading in the M w 8.6 event by about 0.25 bar in the study area. Steacy et al. (2005) pointed out that stress transfer leads to changes in the probability of earthquake occurrence for large earthquakes and aftershocks. Therefore, the phenomenon of increasing SHF or probability of exceedance of earthquake occurrence in this study is then associated with previous studies on the CFS. We interpreted that the increasing SHF obtained in this study is most likely due to the increase in the region's stress load, which was well documented by the GPS data. As shown in Figure 6, increasing the CFS by about 0.25 bar (Pollitz et al., 2006) could increase the PGA level by nearly 0.1. The result may be very beneficial for earthquake mitigation and modeling efforts for probabilistic seismic hazard study and future analysis, especially in the Sumatra area.

CONCLUSION
The SHF analysis off the west coast of northern Sumatra is challenging since there is a lack of information regarding geology, geodesy, and seismological records of seismicity. The result shows that SHF based on the horizontal crustal strain data model tends to be higher than the SHF based on the earthquake catalog and the SHF result based on geological rate. We interpreted that the increasing SHF is most likely due to the increase in the region's stress load. It is suggested that it was the increased stress loading on December 26 of the 2004 event. This is confirmed by the results of a previous study on CFS, in which it was pointed out that the M w 9.1 event caused an increased CFS on the M w 8.6 event by about 0.25 bar (Pollitz et al., 2006). This study's critical finding is that increasing static stress could increase the probability of exceedance of the PGA level in the SHF compared. The increase of the CFS by about 0.25 bar (Pollitz et al., 2006) could increase the PE of the PGA level by nearly 0.1. Owing to the wealth of GPS station data available on northern Sumatra, which recorded the pre-and co-seismic events in 2005 around the study area, the seismic hazard function around off the west coast of northern Sumatra could be well constrained and analyzed. Evidence from the Nias-Simeulue area is an example that clearly shows the effect of increasing static stress on the SHF. It is in line with what was pointed out by Petersen et al. (2004) regarding the possibility of a large increase in the potential for earthquake hazards between the Nias and Simeulue Islands. This increased hazard is based on an M w 8.6 earthquake that occurred in 2005. The location of that event is very well constrained using GPS data distributed in the study area. Thus, this model represents the most updated seismic hazard function in the region. The results obtained in this studymay be very beneficial for earthquake mitigation and modeling efforts for probabilistic seismic hazard study and future analysis.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/supplementary material, further inquiries can be directed to the corresponding author/s.