Site Coefficient of Short Fa and Long period Fv Maps Constructed from the Probabilistic Seismic Hazard Analysis in Yogyakarta Special Province

Among the devastated earthquakes that have occurred in the Yogyakarta Special Province (YSP) were earthquakes that occur in 1867, 1943 and 2006 [1,2]. Macro site coefficient is already available in [3] however, more detailed site coefficient or presented on a micro scale is not yet available. The disaster mitigation program is needed and one of which is the availability of the microseismic site coefficient maps for short period Fa and long period Fv particularly for building design. The probabilistic seismic hazard analysis (PSHA) has been carried out. The subduction and the shallow crustal earthquakes within 500km radius from the city of Yogyakarta have been used. In contrast to [4], PSHA was carried out by using 10% probability of the earthquake being exceeded within 50 years. Result of this study indicates that YSP is dominated by hard SD and medium soil categories. The medium soil is mainly located at along the Opak and Progo river valleys and along the southern coast of Bantul and Kulonprogo districts. The short period site coefficient Fa is strongly influenced by the Opak fault as an earthquake source, while the site coefficient for long period Fv is more influenced by the Megathrust earthquake source. For PGA < 1.0g, the seismic coefficient Fa for short period increases with increasing Vs30, meanwhile for long period of seismic coefficient Fv decreases with increasing Vs30.


Introduction
There have been a number of earthquakes in Yogyakarta Special Province (YSP) that caused great damage and loss of wealth. According to some sources, in the region of Yogyakarta has experienced many occurrences of earthquakes including in 1840, 1867, 1875, 1943, 1947, 1981, 1992, 2001, 2004 and 2006. Based on the data, the Yogyakarta earthquake which caused the most severe impacts was the earthquake of 1867, 1943 and 2006. All of them were subduction earthquake, except the Mw6.3 of May 27th, 2006 Yogyakarta earthquake which was categorized as a shallow crustal earthquake.
The Yogyakarta earthquake on June 10, 1867 [1] caused considerable human damage and casualties, especially in Yogyakarta, Surakarta and south of central Java. The earthquake intensity measured by the Modified Mercalli Intensity Imm scale reached Imm = VIII to IX. The death toll in Yogyakarta is reported to be more than 500 people, more than 440 stone houses completely collapsed. Meanwhile, the death toll in the city of Surakarta were more than 230 people and there has also been a large damage to buildings. The tremor of the earthquake was reportedly felt to the islands of Bali and West Java.
Another major earthquake that occurred in Yogyakarta was on 23 July 1943 with coordinates of the epicenter of 8.6S, 109.9E The earthquake had a magnitude of M7.1, with a depth of 60 km and the coordinates of the epicenter were 8.6S and 109.9E [5]. The earthquake had caused as many as 213 fatalities and more than 12,000 houses were damaged starting from Garut, south coast of Central Java, Yogyakarta, and Surakarta. It was also reported that in Bantul there were 31 fatalities, 654 injured and more than 2600 houses collapsed completely [1].
The 2006 Yogyakarta earthquake with a magnitude of M6.4 actually still included moderate earthquakes and did not include a strong earthquake. According to USGS, Yogyakarta earthquake epicenter coordinates are 7.96S, 110.45E, with earthquake epicenter on land or also called inland earthquake The epicenter is only about 20 km from the city of Yogyakarta with an earthquake depth of only 12.5 km [6]. According to statistical data, the average population density of Yogyakarta is 585.6 people / km2, so this is in the category of very dense population. With such conditions, the death toll caused by Yogyakarta earthquake is more than 5500 people with more than 150000 houses totally collapsed [2].
With the earthquake events as mentioned above, the disaster mitigation program particularly the earthquake natural disaster in YSP is very essential. Many earthquake disaster mitigation programs can be carried out, one of which is the provision of seismic hazard maps through Probabilistic Seismic Hazard Analysis (PSHA). By using PSHA the important parameters in earthquake resistant building design namely the site coefficients at the short period (0.2s) Fa and at the long period (1.0s) Fv can be determined.
The seismic site coefficient of the soil layer has been discussed by [7]. The frequency domain (FD) equivalent linear (EQL) and time domain (TD) non-linear (NL) analysis is the common approach used for performing 1-D seismic ground response analysis. The study was carried out by assuming soil layer thicknesses were 30, 100, 300, 500 and 1000m respectively and its corresponding soil modulus reduction curve and damping curve were adopted from [8]. The seismic soil factor is computed by comparing the spectrum acceleration at the ground surface and the reference/bed-rock ground motion. Results of the study indicated that the seismic site factor/coefficients computed from frequency domain (FD) equivalent linear (EQL) quite different from the time domain (TD) non-linear (NL) site response analysis. In addition, the value of the site factor is also significantly affected by the thicknesses of the soil layer.
The estimation method in determining the seismic site coefficient has been carried out by [9]. It can be categorized as an estimation method because the exciting force used in bedrock is a synthetically generated earthquake time theory through stochastic finite-fault method. The site response analysis is carried by propagating the estimated bed-rock acceleration to the ground surface by using the SHAKE-91 package program with an equivalent non-linear approach. Site coefficient for both short period Fa and long-period Fv was obtained by comparing the ground surface to bedrock response. In addition, uniform hazard spectrum (UHS) in bedrock has also been obtained through the PSHA process. Another result is smooth surface response spectrum (SSRS) obtained through a product between UHS and Fs site coefficient Development of land class maps and site coefficients in the Semarang area by applying one of the estimated methods has also been carried out by [10]. Two PSHA models have been carried out by applying two different shear wave velocity VS values. The first model or M-1 is PSHA which uses shear wave velocity in bedrock whose value is Vs = 760m/s. The second model or M-2 is PSHA on the ground surface based on the value of the nearsurface shear wave velocity at the uppermost 30m soil layer or VS30. The mechanism of the earthquake source used was a the subduction earthquakes that occurred at the Indian Ocean south of Java island and shallow crustal earthquakes due to the fault activity at a radius 500 km from the city of Semarang. The results showed that the Fa and Fv site coefficient obtained from the study was relatively close to the Fa and Fv written in the Code. The main aim of this study is to develop the site coefficient maps at Yogyakarta Special Province (YSP) by applying the combination of the method that has been used by [9] and [10]. In the meantime, the shear wave velocity data in YSP is not yet available therefore the VS30 data published by USGS has been used. Probabilistic seismic hazard analysis (PSHA) based on the total probability theorem with seismic data from USGS has been used where the analysis was based on 3-D earthquake source.

The Geology of the Yogyakarta Special Province (YSP)
As shown in Figure 2), the Yogyakarta Special Province (YSP) consists of 4-districts and 1 municipality. In the north is Sleman regency, in the south-east is Gunungkidul district, in the south is Bantul regency and in the west is Kulonprogo district. Also seen in the picture is that Mount Merapi is located at the northern end of Sleman district. Meanwhile, the southern part of YSP is the Indian Ocean where the Australian Plate moves towards the north and subduction occur approximately 250 km south of the YSP coast. On the other hand, the island of Java is located on a Eurasian plate that moves southward so that it collides with the Australian plate on the boundary plate. With the convergent movement of tectonic plates to each other, depression areas have occurred which will further affect the geological conditions at YSP.
According to [11], the formation of the rocks and soils in YSP were initiated during the Tertiary and Quarternary periods. As shown in the Fig.2) the western part of Kulonprogo and Gunungkidul districts are mountain ranges where the rocks were formed during the Tertiary period. In the northern part, i.e the Sleman regency (the closest to the slopes of the active volcano Mount Merapi), the soil profile consists mainly of coarse sand and gravel,  [12]. The soil in these areas consists of relatively thick and young soil deposits, which have a high potential seismic vulnerability index [13].

Earthquake source mechanisme of YSP
It was stated earlier that the java island including the YSP located in the Eurasian Plate moved southward colliding with the Australian plate which is moving northward. Therefore there is a subduction of approximately 250 km to the south of the island of Java where the Australian plate subducts under the Eurasian plate. In the subduction area, there is an interface-slip between two tectonic plates which is the source of the earthquake both past and future. The earthquake source mechanism in the subduction area then is very commonly referred to as Megathrust and Benioff zone. The earthquakes that occurred in Megathrust and Benioff zone became one of the sources of the earthquake in YSP.
The earthquakes that occurred in the Megathrust zone in the southern Java island were likely to be more than Mw8.0. However, because of its distance is more than a couple hundred kilometers, then the effect on ground acceleration in Yogyakarta is relatively  small. The earthquakes that occur on land or inland earthquakes are earthquakes that occur due to fault activity. Although the magnitude is not as big as the Megathrust earthquake, because of its distance is close and the hypocenter is not deep, the effect on ground acceleration is likely to be greater than subduction earthquake. Therefore, in addition to subduction earthquakes, shallow earthquake or shallow crustal earthquakes which are caused by fault activities will be the source of the earthquake at YSP.

Probabilistic Seismic Hazards Analysis
In the probabilistic seismic hazard analysis (PSHA) the probability of an earthquake event with parameter X that exceeds an x value for an earthquake with magnitude m can be written to [14,15], (1) when the earthquake distance R is considered, then the probability of hazard parameter X exceeds an x value for a given magnitude M and distance R can be expressed in an integral form as, M, m is earthquake magnitude, R, r is the earthquake distance, fM and fR respectively are the probability density function of magnitude and distance.
Meanwhile, the annual rate of exceedance for earthquake X parameters for a particular earthquake source can be determined through, (3) where N is the average number of earthquake events per year for earthquake magnitude ranging from mo to mmax obtained from one of the earthquake sources in the field.
For some potential earthquake sources that will affect seismic hazard at the site being reviewed, the annual rate of exceedance can be determined in the same way. If the number of potential earthquake sources is i = 1 ... n, then the total annual rate of exceedance lX(x) is the contribution of all earthquake sources that can be written into, for a 3-dimensional earthquake source, the probability of the occurrence of parameter X that exceeds the x value must also take into account the influence of attenuation. Therefore, the effect of earthquake magnitude, distance and attenuation must be combined, so that Eq.4) becomes, the annual rate of exceedance as presented in Eq.3) is calculated based on the pseudospectral acceleration (PSA) for a given T vibration period, so that a hazard curve can be made for each desired T period. Based on the hazard curve then uniform hazard spectrum can be made for each point under consideration. Thus bed-rock, short period and long period hazard map can be made.

The Site Coefficient, Fa and Fv
The seismic soil site coefficient is defined as the ratio of the spectral acceleration at the site with its reference rock outcrop at the same period [16]. The empirical site coefficient decreases with increasing the level of ground acceleration. The empirical value of site coefficient is commonly presented in short periods (T = 0.2 s) and long periods (T = 1.0s) which is generally abbreviated as Fa and Fv. There are several concepts or approaches to determining seismic site coefficient. As stated by [7], the site coefficient is influenced by the thickness of sedimentary soils. According to [17,18], the seismic site coefficient at short period Fa and long period Fv can be determined by, where Rsoil and Rrock are the hypocenter distance of soil and rock respectively, RSsoil and RS rock are the response spectrum at ground surface and reference point/bed-rock and T is period.
Meanwhile according to [19], the site coefficient at short period Fa and long period Fv can be determined by the following expressions, It can be identified from Eqs.8) and 9) that aside from being affected by the thickness of sedimentary soil layer, the short and long period site coefficient is affected by the shear wave velocity at uppermost 30m or VS30. The approximate value of Fa and Fv as a function of VS30 have been proposed by some researchers.

Earthquake Data
As stated earlier that the seismic source mechanism identified for the YSP is subduction and shallow crustal earthquakes. Subduction earthquakes were separated into 2 zones, namely earthquakes in the Megathrust and earthquakes in the Benioff zone. Meanwhile, shallow crustal earthquakes are earthquakes that occur around the fault. Earthquake data from both subduction and shallow crustal from 1963 to 2017 within a radius of 500 km from the city of Yogyakarta are considered. The earthquake data taken were earthquakes with magnitudes greater than M5.0 with a depth of less than 700 km. The earthquake data was downloaded from the USGS catalog.
Earthquake data collected is still general, meaning the earthquake data is still mixed including foreshock, main-shock and aftershock earthquakes. The earthquake data needs to be processed statistically, earthquake foreshock and aftershock are separated from the main shock, furthermore, only main-shock earthquakes are used for the analysis process. The main-shock data was then transformed into an earthquake with the same magnitude, namely the moment magnitude Mw.

Earthquake source identification
According to [1] the Cimandiri fault located at West Java as presented in Fig.2) is a strikeslip fault. The Cimandiri fault with a strike of N70o-80o, extends from the Gulf of Pelabuhan Ratu and passes through the Sukabumi, Cianjur and Bandung area. The Cimandiri fault is about 100 km length and undergoes slip-rate about 0.5-1.7 mm/year by [20], 4mm/year by [4]. Meanwhile, the Bumiayu fault that is located in the western part of Central Java is an extension of the Baribis fault. According to [21], the fault of Bumiayu is right lateral strike-slip with dip angle 90o and slip rate of 2.0 mm/year, and has caused an earthquake in 1995.
In Central Java, the Lasem fault is an active fault starting from the eastern part of the city of Semarang towards the east through the northern mountains of Purwodadi and continues up to Lasem [10]. Meanwhile, the Pati fault is a branch of the Lasem fault, continuing towards to the east at southern of Kudus and Pati town until near to Rembang. Data and parameter of Lasem and Pati faults are presented in Table 1.  Fig. 3. The earthquake source mechanism of the YSP.
In YSP, the Opak faults were not well known before May 27th, 2006 Yogyakarta earthquake. According to [22] the 27th may 2006 Yogyakarta earthquake was caused by the activity of oblique Opak fault with strike-slip of 0.8m, dip-slip of -0.26m and strike N 48o to E. The slip-rate of Opak fault is about 2.4 mm/year with the maximum historical earthquake magnitude of Mw6.8 [4,5], According to [23] the Opak fault has a dip angle of 87 o ; rake angle of 177 o (right lateral strike-slip) and the upper edge of the fault (ZTOR) at 6 km. The Opak, Lasem, Pati, and Bumiayu faults are the active faults that are taken into account in this analysis.
Another earthquake sources mechanism that must be considered in YSP is source caused by the subduction of Australian with respect to the Eurasian plate. As shown in Fig.  3) and Table 2, within a 500km radius of the city of Yogyakarta, the subduction earthquake source was divided into 3 zones i.e. zone-1, zone-2, and zone-3. Each zone consists of Megathrust and Benoff earthquake sources. According to [4,21], subduction activity in the south of the island of Java has a potential to cause a maximum earthquake of Mw8.1.

The Ground Motion Prediction Equation (GMPE)
In general, GMPE an earthquake parameter is stated in an equation that contains 3 main components, namely the source, transmission path and site effects [24]. In the beginning, the GMPE was only presented in terms of peak ground acceleration (PGA). In the next development, the GMPE is presented in the period vibration T function to produce spectral displacement (SD) and pseudo-spectral acceleration (PSA) or pseudo spectral velocity (PSV) as a simplification. The researchers suggest that the selection of GMPE must be done carefully and appropriately. The general and detailed criteria for selecting and adjusting the GMPE models have been presented briefly by [25,26]. Among the general criteria is to consider global earthquake models and exclusion criteria. According to [25,27,28], the exclusion criteria in the selection of the GMPE it should be considered in the PSHA and the criteria have been summarized by [29]. As stated earlier that the source of the earthquake in YSP is mainly generated from the subduction and shallow crustal mechanisms. Therefore, in the PSHA several GMPEs were chosen which were not included in the exclusion criteria as presented by [ ] respectively for the subduction and shallow crustal earthquake sources Considering that GMPE derived based on earthquake data in Indonesia is not available, then the GMPE from other countries is used where the earthquake source mechanism is similar to in the YSP. The GPME used for modeling the parameter of shallow crustal earthquake in YSP are: a) Boore, Joyner and

The Logic Tree
Researchers have identified that when GMPE is formulated, there are many aspects of uncertainty that can reduce the accuracy of GMPE. Based on many studies that have been carried out, in general, the researchers agree that there are intrinsic uncertainties namely sand aleatory uncertainty. [30][31][32].The epistemic uncertainty can be caused by an incomplete knowledge which can consist of several things i.e. inexact model selection, error in statistics, error in measurements and errors in the data base. In probabilistic seismic hazard analysis (PSHA) the presence of epistemic uncertainty can be treated by applying the logic tree model [33]. Meanwhile, aleatory uncertainty in another side may be caused by the natural variability of the soil site, the soil non-linearity or site amplification [31] usually consists of intra and inter-events [34] and presented by standard deviation or sigma [25].
In the logic tree, there are several GMPE that can be used, each of which has a different degree of uncertainty. Because of the degree of uncertainty, it eventually leads to the use of the weighting factor for each GMPE used. To increase objectivity and ranking, each GMPE candidate is assessed through expert judgment and tested through procedures such as presented by [35]. Finally, a weighting factor is given to GMPE candidates based on the result of ranking and expert judgment that has been decided. Based on all important aspects as mentioned, the logic tree for shallow crustal earthquake source in the PSHA at the YSP is as presented in Figure 4.

Site Classes
It is common that natural soils can consist of several types such as cohesive soils, non-cohesive soils or rocky soils. Each type of soil has its own static and dynamic properties which will then affect the soil behavior due to an applied load. Both static and dynamic, each soil property will have its own parameters both physical and non-physical parameters. Under the dynamic loads due to the earthquake, for example, each type of soil will have a basic response behavior that will depend on the loading intensity, type and internal properties of the soil. In earthquake engineering, soil types that are still on a micro scale are then extended to a macro classification which later became known as site classes. Site classes in YSP and its surroundings that are presented in the form of maps are as shown in Fig.5. Site classes are based on VS30 data downloaded from the USGS database grouped according to [36] which is based on the near-surface shear wave velocity VS30. According to [36]

Site Coefficients or Fa and Fv Maps
It is generally known based on earthquake recordings in the past and through site response analysis that soil amplification has occurred in thick sedimentary soils due to an earthquake. For earthquake resistant design of the building, soil amplification is needed at the building site, but site response analysis cannot always be carried out for several reasons. Therefore in each Building Code, the value of site coefficients Fa and Fv has been given through macro analysis as an approximation of soil amplification.
The site coefficient for short period T = 0.20s or Fa for pseudo spectral acceleration (PSA) less than 1.0g in the YSP and surrounding is presented in Fig.6). The site coefficient was derived based on PSHA with a 2% annual rate of exceedance over a period of 50 years. In the figure, it appears that site coefficient for short period is divided into 4 groups. It appears in the picture that the lowest value of Fa = 1-2 is located along the Opak fault. This is logically caused by the large PSA at base rock near the fault, while the PSA distribution on the ground surface is not the same as in the base rock. Therefore the Fa value near the fault becomes the smallest and enlarges at a farther distance from the fault. It is also clear from the picture that most YSP regions have Fa = 1-2, especially in areas with site class D. Meanwhile, site coefficient for long periods T = 1.0 s or Fv presented in the form of maps is as shown in Fig.7). If Fig. 7 and Fig. 8 are compared, it shows that the distribution of Fa and Fv values is very different. It is clear from the figure there is the only 2-main group that most of the YSP region, especially the northern region, has a site coefficient of Fv = 2-3, while the southern region has a larger site coefficient, Fv = 3-4. These results may be due to the distance effect of the Megathrust earthquake mechanism, because the high value of Fv along the southern coast of YSP is closest to the Magathrust earthquake source.

Site Coefficient as Function of Shear Wave Velocity, VS30
As stated earlier that the value of the site coefficients Fa and Fv is influenced by several variables including the PSA value, Vs30 and the thickness of sedimentary soil [7]. The site coefficients are also can be determined by the ratio of the spectral area as presented in Eqs. 6) and 7). Meanwhile, the site coefficients as presented in Eqs.8) and 9) are a function of shear wave velocity. The calculated site coefficient covering 416 data point in the YSP for PSA < 1.0g by applying the method as presented by [10] is depicted in Fig. 8).
It can be seen from the picture that the site coefficient for long periods of Fv (T = 1.0s) was greater than site coefficient for short periods of Fa (T = 0.2s). This is identical to the site amplification where amplification in ground motion for low frequencies will be higher than site amplification due to a high-frequency ground motion. This happens because of several reasons: 1) not all GMPE used have considered Vs30 as a variable; 2) site amplification besides being affected by Vs30 is also still affected by the thickness of the sedimentary soil and level of applied ground motion.
The site coefficients both Fa and Fv as a function of Vs30 as shown in Fig. 8) appear to be spread widely, but in general, they have a similar pattern. It also appears in the figure that Fa increases with increasing Vs30, while Fv decreases with increasing Vs30. Although the site coefficients distribution pattern in this study is not the same as that presented by [37], in general, the distribution of site coefficients is spread out relative to the Vs30. In the future, broader and deeper studies are needed to estimate the relationship between site coefficients and near-surface shear wave velocity Vs30.

Conclusions
Studies on seismic site coefficient for short period Fa and long periods of Fv in YSP have been conducted. The values of Fa and Fv are estimated through the ratio of PGA on the ground surface relative to PGA in base-rock. The PGA at the ground surface and bedrock is calculated based on the PSHA with a 2% annual probability rate of exceedance for 50 years. The Fa and Fv are presented in the form of a maps and graph with respect to the Vs30. Based on the study, the results can be summarized as follows.
1. According to the USGS database, soil category SD or medium-soil is located along the Opak and Progo river valley, along with the coast of Bantul and Kulonprogo districts and in the middle part of Gunungkidul district. Meanwhile, the hard soil type or SB is located on the slope of Merapi, the western part of Kulonprogo district and along the edge of Gunungkidul district.
2. The site coefficient at a short period or Fa is affected by the sedimentary soil at Opak river valley or strongly affected by the activity of Opak fault as a seismic source. For PSA < 1.0g, the farther distance parallel with Opak faults, the site coefficient for the short period Fa will be even greater. Meanwhile, site coefficient for long period Fv seems to be strongly influenced by earthquake activity in megathrust. This happened because a large Fv value occurred along the southern coast of YSP parallel to the Australian plate boundary, 3. The relationship between site coefficients Fa and Fv with Vs30 seems to be scattered. It was found that the seismic coefficient Fa for short period increases with increasing Vs30, meanwhile for long period of seismic coefficient Fv decreases with increasing Vs30