Development of synthetic ground motion at a specific site in Yogyakarta town, Indonesia utilizing the PSHA Method

: This paper presents the development of synthetic ground motion at specific sites in Yogyakarta town. In the 2019 Indonesian Seismic Code [1] provides an alternative method in the analysis of building structures by applying the dynamic time history analysis. At least 11-pairs of earthquake recordings must be used in the analysis. Synthetic ground motion utilizing the Method of Probability Seismic Hazard Analysis (PSHA) was used in this study. A selected site in Yogyakarta town was chosen as a pilot study considering that there were many fatalities and building damage caused by the 2006 Yogyakarta earthquake. The Uniform Hazard Spectra (UHS) based on the shallow crustal earthquake source is higher than the Megathrust. The risk targeted spectrum demand MCEr has been considered, which on average 12.3% greater than the UHS. The synthetic ground motions (SGM) are accordingly based on the shallow crustal earthquakes. The dominant magnitude and distance are M D = 6.5 and R D = 14.5 km. They show that the contribution of the Opak River fault to the hazard in Yogyakarta town is very dominant because the distance is very close. Based on the obtained M D and R D , spectral matching, and testing significant duration D 595 , the 12-synthetic ground motions were successfully developed.


Introduction
Indonesia is one of the earthquakes (EQ) prone countries in the world, including in the Yogyakarta Special Province (YSP).As stated by Pawirodikromo [2], there is no significant strong-motion record network in Indonesia, so the availability of ground motion time series is still a big problem.The problem was related to structural analysis after the 2019 Indonesian Seismic Loading Code [1] was published.Indonesia needs to build a network of reliable seismic records, such as being installed in Turkey [3], China [4], or Taiwan [5,6].
Besides the spectrum response, the structural analysis in the 2019 Indonesian Seismic Loading Code [1] allows using the time history analysis methods (based on both elastic and inelastic responses).If the Performance-Based Seismic Design (PBSD) is used, then a large number of ground-motion time series are required [7], similar as needed in the 2019 Indonesian Seismic Loading Code.In the 2019 Indonesian seismic Code, at least 11-pairs of earthquake time history must be used, both from Megathrust, Benioff, or shallow Crustal source mechanism.The reliable availability of earthquake time series for site analysis or structural response is one of the biggest and complex problems in earthquake engineering.It happens because many earthquake source variables included, and many site uncertainties occur [8].
According to Bulajic and Manic [9], Rezaeian, and Kiureghian [10], the availability of SGM can be created by using the stochastic, deterministic, and probabilistic methods.In this paper, the SGM was created by using the Probabilistic Seismic Hazard Analysis (PSHA).One of the PSHA results is the availability of seismic hazard maps, in particular, are under consideration.
The need for seismic hazard maps has long been initiated by experts, especially after the introduction of the engineering seismic risk analysis concept by Cornell [11] for line and area source, particularly in the USA.Meanwhile, Atkinson [12,13] said that at almost the same time, the probabilistic seismic hazard was developed in Canada.Recently, 3-D earthquake sources or fault sources have become standard models used in computation routines in PSHA, including seismic hazard maps that have been published in the 2019 Indonesian Seismic Loading Code.
As stated earlier, 11-pairs earthquake time series as required can be generated through two approaches: 1) through the matching of the spectral design from the Code, and 2) through shear wave propagation in soil media of the SGM in bedrock.In term of macro zonation purposes, the first approach has a weakness, since it can only be used in certain regions spatially and cannot be used for specific-site analysis.For specific site analysis, in contrast, earthquake ground motion is needed from PSHA in which the SGM on the bedrock level will be obtained.Besides, the second approach requires a more significant effort than the first approach.This approach happens because the latter method requires double spectral matching.The first matching in the second approach is spectral matching for each SGM at bedrock level to UHS as obtained from PSHA.Meanwhile, the second matching is matching of the earthquake time series on the ground surface (as a result of the wave propagation process of the SGM at bedrock) against the mean spectrum of earthquake time series at ground level.Considering that the second stage is very long and time-consuming; therefore, in this paper, the discussion is limited to the first phase, i.e., the earthquake timeseries matching on the bedrock level.Based on those as mentioned above, then the objective of this study was to obtain at least 11 SGMs in a bedrock level that were most likely to occur at one selected point in Yogyakarta city.

A Previous Study of the SGM
An example of the survey of SGM in Yogyakarta was carried out by Sunardi [14].There are three critical elements in generating a matched spectrum that is targeted range, original/artificial ground motion, and matching process.In this study, the targeted spectra were created through the PSHA standard.The recorded earthquake ground motion was selected from the USGS catalog, while the matching process was carried out by using the standard Fast Fourier and Inverse Transform method.Deaggregation in the period T = 0.2 s and long periods T = 1.0 s were used.Although the generated SGM based on 3-dimension of earthquake sources were used, the generated SGM is limited only to one motion in each earthquake source and does not discuss the seismic intensity measures.Thus, the result of this study does not meet the requirements of 11-earthquakes as required by [11].
A comparison between the real earthquake and SGM has been presented by Parajudi and Shrestha [15].In this study, two target spectrums were used based on the deterministic method, i.e., by using Young et al. (1997) ground motion prediction equation (GMPE); and the Barpak Nepal earthquake.Meanwhile, the artificial earthquake is generated using the stochastic method.The results of the study showed that the generated SGM was quite close to the Barpak real earthquake records.This study is also limited to the discussion of 11 earthquake records.
An intensive demonstration of the compatibility between the response spectrum from artificial ground motion (generated by the stochastic method) to the target spectrum was carried out by Huang and Wang [7].Besides matching the spectrums, the developed method also used for matching the energy content (energycompatible) through Arias Intensity seismic measure.The whole process is called energy compatible and spectrum compatible (ECSC) and carried out through the iteration scheme.In this study, 50 simulated ground motions were generated and compared to 50 real recorded earthquakes (PEER-NGA).The 50 simulated earthquakes and its counterpart of 50 recorded earthquakes were utilized in the inelastic dynamic analysis of 12-story reinforced concrete frame building.The results of the study demonstrate the structural response due to ground motions generated by ECSC and PEER-NGA records are very consistent.This method has successfully created and simulated its compatibility between simulated and real motions; however, it does not generate the target spectrum for design purposes.

Earthquake Sources of the YSP
The YSP is located in the south of Central Java Province and is directly adjacent to the Indian Ocean, as shown in Fig. 1.Tectonically, YSP or Java Island is located in the southern part of the Eurasian Plate which moves towards the south [16], and opposite collides in long subduction plate boundary with Australian Plate which is moving to northward with the average movement of ± 5.5 -5.7 cm/year [17,18].With such conditions, the Megathrust earthquakes in the plate boundary in the south of Java cannot be avoided in which an example of the crosssection of earthquake events is as presented by [19] or as depicted in Fig. 2.    In addition to the Megathrust earthquake, in the YSP area was also shaken by the M7.2 (coordinate of 8.88 S and 110.65 E) and the M8.1 Benioff earthquake in 1943 with coordinates of 8.6 S and 109.9E with a depth of 90 km [20].As illustrated in Fig. 2, the plate boundary in the south of the island of Java is located approximately 325 km south of the city of Yogyakarta.By considering the distance, thus, the effect of the Megathrust earthquakes on the city of Yogyakarta is estimated to be relatively small.In addition to the impacts of Megathrust earthquakes and Benioff earthquakes, the YSP area is also affected by shallow earthquakes both due to the activity of collision plates and due to active fault activity called the Opak river fault.According to Pawirodikromo [2], the Opak riverfault, which is about 54 km long, is located approximately only 10 km south-east of Yogyakarta city.

2006, Epicenter
The study of the mechanism of the Opak fault movement has been carried out by many researchers since the occurrence of May 27, 2006, the Yogyakarta earthquake, where the results differ from each other.Opak faults are modeled by Thant et al. [21] moving with a dip-slip mechanism with SW-NE orientation.Meanwhile, the Opak fault movement is modeled by Natawidjaja [22] with left-lateral strike-slip, and this mechanism confirms with the results study, as reported by Abidin et al. [23].The results of their study show that the Opak fault is N48 o E oriented with 89 o dip slip.Meanwhile, Tsuji et al. [24] interpreted that 2006.Yogyakarta earthquake as a result of the oblique mechanism movement on the Opak river fault and left lateral strike-slip on the fault rupture that is estimated to occur at 10 km east of the Opak River.

Annual Rate of Exceedance and Deaggregation
The discussion and development of PSHA have received much attention by researchers and has been presented on many occasions since the application of seismic risk in engineering project introduced by Cornell [11].Much overall understanding, refinement of the method, and its application in Seismic Codes have been widely used in many countries.A prescribed hazard essentially needs to be integrated for all possible earthquake magnitude, distance, and earthquake source through the total probability theorem through the formula [12] to determine the probability of exceeding a specified ground motion amplitude, where  (M > m min ) is the annual earthquake rate of occurrence greater than minimum magnitude m min , P(IM >x|m,r) is the probability of intensity measure IM greater than x caused by variables earthquake magnitude m and distance r, F M and F R respectively are cumulative probability function of m and r.
When the SA values in each period T have been obtained, the results are then presented in terms of UHS.Each SA value in the UHS graph is the result of contributing to all possible types of sources, earthquake distance, and magnitudes.It is essential, therefore, to know the dominant magnitude (M D and distance (R D ) that make a dominant contribution to each value of SA and it respectively presented as [25,26]: M i,D is the earthquake magnitude that contributes most dominantly to hazard at one point due to all possible magnitudes and distances caused by i th earthquake source,  i,m,r (x) is the annual rate of exceedance due to the i th earthquake source that has taken into account all possible distances and earthquake magnitudes,  m,r (x) is an annual rate of exceedance due to all earthquake sources and has taken into account all possible distances and earthquake magnitudes.Meanwhile, R i,D is the earthquake distance that contributes most dominantly to hazard at one point due to all possible magnitudes and distances caused by i th earthquake source,  i,m,r (x) is the annual rate of exceedance due to the i th earthquake source that has taken into account all possible distances and earthquake magnitudes,  m,r (x) is an annual rate of exceedance due to all earthquake sources and has taken into account all possible distances and earthquake magnitudes.

Synthetic Ground Motion (SGM)
The first output taken from the PSHA results is the relationship between the spectral acceleration (SA) and the annual rate of exceedance that is displayed graphically in the form of a hazard curve.For example, a yearly rate of exceedance of 0.2% is taken (10% probability of exceedance for 50 years), then from the hazard curve will produce a relationship between periods T and spectral acceleration (SA), which is then called the UHS and functions as a target spectrum of the site.On the other hand, the value of the spectral acceleration in the UHS has not been able to provide information on the magnitude and distance that contribute most dominantly to spectral acceleration.For this reason, it is necessary to do a re-breakdown or deaggregation to obtain the dominant magnitude and distance (M D and R D ).Deaggregation can be done for a combination of all earthquake source mechanisms or separated for each type of earthquake source.
In practice, a local earthquake record that matches the obtained M D and R D is not necessarily available.For this reason, it is necessary to look for records in the earthquake Catalog that match or are close to M D and R D .
Based on the earthquake record, then the actual spectral acceleration can be generated.Thus there are two spectra, namely the target and the real spectral acceleration, which is not necessarily matching.Therefore, the next step is to do spectral matching, which will eventually result in matched synthetic ground motions on the bedrock.

The D595 testing of SGM
This study is related to the availability of SGM to one point in the city of Yogyakarta due to various earthquake source mechanisms.However, because the river fault Opak is the closest earthquake source, then the shallow crustal becomes the most dominant earthquake source.The study then leads to the site-specific response, such as carried out and its benefits by Bradley [27].Since the selected point/site and earthquake source are already more specific, then the possibility of a dominant earthquake that will occur and recorded in the site is not expected very random but is likely to be relatively close to one another.The proximity value criteria proposed in this study is the earthquake significant duration D 595 as proposed by Kempton and Stewart [28] and Lee et al. [29] and expressed in the equation,  (4) in which the D 595 is the earthquake significant duration M w is moment magnitude, M r is reference magnitude (M r = 6), V SR is shear wave velocity (m/s), b 1 , b 2 , c 1 , c 2 are coefficients equal to 2.79, 0.82, 0.15 and 1.91.
The significant duration D 595 is commonly used and defined as the length of interval time of which the dissipated energy within 5 -95 % of the total energy of an earthquake ground acceleration [28,29].Thus, the obtained ground motions must have a significant duration within a specific range, according to Eq. 4.

GMPE in the Next Generation Attenuation (NGA)
version is a prediction of ground acceleration at the average or geometric mean level.Therefore UHS obtained is also only in the average or geometric mean spectral demand.Calculating the maximum value in the directivity effect is necessary with regard to computing the value of maximum spectral demand so that it will reach the level of Maximum Credible Earthquake (MCE).Between geometric mean and maximum spectral demand is connected by a value called the directivity factor, which is the function of the period T [30].The directivity factor for T = 0s can be taken D f = 1, for the period T = 0.2s the value of D f = 1.1 and for the period T = 1.0s, the value of the directivity factor D f = 1.3.From T = 0s to T = 0.2s and from T = 0.2s to T = 1.0s, the D f value can be calculated linearly.For T > 1.0s, the directivity factor is constant, D f = 1.3.
Apart from the value of the directivity factor, it is still necessary to know that at the MCE level, it cannot guarantee the occurrence of uniform building collapse during 50 years [31].Therefore, in order to uniform risk to occur, one more factor needs to be taken into account, namely, the risk targeted coefficient (R f ), so that the ground motions will reach the Risk Targeted Maximum Credible earthquake MCEr.Based on the Indonesian Seismic Code, the risk targeted coefficient has been presented in the form of a map that is a C RS map for a short period T = 0.2 s and C R1 map for the long period T = 1.0s.

Method of Investigation
The reason why the location of the study was conducted at YSP was that the 2006 Yogyakarta earthquake had caused many fatalities and property loss.The research that has been carried out using the method or steps as presented below.

Collection and Processing of Earthquake Data
The data used in this study are earthquake selected from the USGS Catalog starting from 1963 -2016 covering an area with a radius of 500 km from the city of Yogyakarta, as presented in Fig. 3.The data is then processed by separating main shocks and aftershocks and converted into the same magnitude unit.

Identification and Modelling of Earthquake Sources
In this step, the earthquake sources are modeled, which are generally classified as Megathrust, Benioff, and Shallow crustal earthquake sources.The characteristics of all earthquake sources are as presented in Tables 1 and  2. The next step is like a routine work that is generally needed before the PSHA process.

PSHA Analysis, Target & Actual Spectra and SGM
The result of this step is a target spectra at a specific point in Yogyakarta city as desired.Meanwhile, the M D , R D (Eqs. 2 and 3), and actual spectra will be obtained during the deaggregation process.The SGM can be generated after matching between the actual and the target spectra.The SH (Seismic Hazard) Model software was used in this study.

D595 Testing
At this step, the D 595 testing is carried out until at least 11 SGM are obtained through the trial and error process.The acceptance test criteria are that the D 595 must be fall in the range of significant earthquake duration calculated by Eq. 4. The lower limit is calculated based on the closest, and the upper limit is based on the farthest epicenter distance.

Uniform Hazard Spectrum (UHS)
Fig. 4 presents a comparison of UHS based on the shallow crustal and Megathrust earthquake sources for one site in Yogyakarta city.It appears in the figure that the UHS due to the shallow earthquake is higher or more dominant than caused by the Megathrust earthquake.The results were obtained because the Opak river fault as a shallow earthquake source is very close to the city of Yogyakarta, while the Megathrust earthquake epicenter as graphically presented in Fig. 2 occurs very far from the city of Yogyakarta.Based on these results, then the earthquake hazard in Yogyakarta that will be discussed further is only a hazard caused by the shallow crustal earthquakes.process show that the dominant earthquake magnitude is M D = 6.5, while the dominant distance R D = 14.5 km.The dominant magnitude M D = 6.5 and dominant distance R D that are so short indicates that the influence of the Opak river fault on the earthquake event in the city of Yogyakarta is powerful.Contributions to hazards due to other active faults such as the Cimandiri, Lembang or Baribis faults are very small because its distance is already very far from the city of Yogyakarta Meanwhile, Fig. 6 is the directivity and risk targeted factors, as stated earlier.It can be seen in the picture that these factors or coefficients are functions or change according to period T. In the picture it appears that the directivity factor D f from T = 0s to T = 0.2s the value of D f > 1 and increases linearly, but for T > 1.0s the directivity factor is constant.Meanwhile, for the city of Yogyakarta the Risk factor value is less than 1.0 or R f < 1.0 for the whole period T. If all these factors have been taken into account, from the hazard spectra demand (UHS) to the risk targeted spectra demand (MCRr) on average, it will increase by 12.3%.Meanwhile, Fig. 7 is a comparative representation between UHS and risk targeted spectra demand (MCEr).It can be seen in the picture that in the entire period T, to the MCEr position, the UHS values have increased by an average of 12.3%, as mentioned earlier.By obtaining risk targeted spectra demand (MCEr), the spectral matching process can be carried out.

Spectral Matching
The MCEr target spectra for a specific site in the town of Yogyakarta has been determined and presented in Fig. 7. Considering there were no recorded earthquakes in the Yogyakarta region [2], then the candidates for earthquake recordings from the USGS Catalog were selected.Finally, 12 earthquake recordings were obtained (more than the minimum requirement of 11 earthquake recordings) as presented in Table 3, all of which are shallow crustal earthquake records.The 12earthquake records came from 5 countries, particularly from the USA, and the others were from Japan, Turkey, Greece, and Italy.
The selection of strong-motion records from the USGS Catalog cannot be carried out arbitrarily, but the frequency content of the spectral target must be considered.Given the spectral target is located at the bedrock level and the dominant earthquake distance R D is very close.Accordingly, the spectral target has highfrequency content.The spectral matching can be taken quickly, the strong motion records from the USGS Catalog with high-frequency content must be chosen.It appears in the table that earthquake magnitude varies from M5.7 -M7.2, which is a typical shallow crustal earthquake magnitude.Meanwhile, the total duration of earthquake records also varied from Tt = 39.93s to Tt = 59.84s.All of the 12-earthquakes as obtained have their actual spectrum, which is not necessarily the same as the target spectra.Therefore matching spectrum needs to be done, which results are as presented in Fig. 8. Fig. 8 is a matching spectrum for the selected 12 earthquake records.In Fig. 8, it appears that in general, 12 earthquakes matching well on the MCEr target spectrum, except at the peak of the spectrum.It is complicated for well matching to the overall spectrum period.Meanwhile, PSV for all 12 matching ground motions is presented in Fig. 9.

Significant Duration of D595 Testing
As stated earlier, the earthquake occurring due to Opak fault activity concerning a specific site in Yogyakarta is expected to be not too random.The significant duration or D 595 then to be a criterion and the results of testing of 12-synthetic ground motions are as presented in Fig. 10.
As shown in Fig. 10 or as presented in Table 3 that the D 595 of the 12 SGMs calculated using Eq. 4 has a range value between 13.04-24.63swith an average of 15.92s.These results still meet within the lower limit of t e = 10.53s and the upper limit of 25.95s.The proportion of significant duration of D 595 = 15.92s is shorter than that of D 595 = 25-30s [33] in the 2016 Kaikoura New Zealand earthquake.Besides, the distance is closer; the maximum earthquake magnitude of Yogyakarta due to the Opak fault is only estimated M = 6.8 and smaller than the Kaikoura earthquake was M = 7.8.
It should be noted that the lower limit t e = 10.53s is calculated based on the closest hypocenter distance to the river the Opak fault while the upper limit te = 25.95s is calculated based on the farthest hypocenter distance, each assuming an earthquake depth of 12 km (such as the 2006 Yogyakarta earthquake).With the fulfillment of all the significant duration, then the 12 SGMs, as presented, can be used as the basis for further analysis.

Synthetic Ground Motion (SGM)
After a spectral matching and a significant duration test of D 595 , the 12 proposed SGM that are likely to occur in Yogyakarta city due to Opak river fault activity are presented in Fig. 11.It appears in the figure that by a glance, the chosen SGM is relatively varied in both the PGA and frequency content (which is indicated by its value of A/V ratio).The picture also presented the accumulated seismic energy, which is mathematically presented by the integral of the square ground acceleration over its duration.It appears in Fig. 11 that the red line in the figures shows one-tenth of the accumulated energy contained in earthquake records.It appears that the El Centro 1940, Hector Mine, Kobe-Ak, and Mammoth Lakes earthquakes contain the highest energy content as compared to the others.It is affected by the PGA value and the duration of the strong part portion on the earthquake record, as clearly shown in Fig. 11.
Note that all 12 SGMs are records in the East-West Component (EWC), while its partner in the North-South Component (NSC) will be presented on another occasion.Once again, the 12 SGMs, as presented, are most likely earthquakes generating only from the seismic activity of the Opak river fault identified at one point in the city of Yogyakarta.Even though the location of the specific site, earthquake source, and wave transmission path to the site are definitive, the possibility of future earthquake characteristics may still be different.Besides the total duration Tt and significant duration D 595 as mentioned, the 12 SGMs also has varied maximum ground acceleration as presented in Fig. 12.
The variation of PGA and frequency content (A/V ratio) of 12-SGMs can also be seen visually in Figs. 12 and 13.Based on Fig. 12, it can be seen that PGA varies from 0.183-0.276gwith an average value of 0.232g.Meanwhile, based on Fig. 13 the A/V ratio also varies from 0.771-2.177,with an average of 1.329.Almost all SGMs fall in the high earthquake frequency content [34] because it has an average A/V ratio higher than 1.20 and confirms to the condition that the earthquakes are still in the bedrock level.

Conclusions
The SGM studies at the bedrock level, which might occur in the city of Yogyakarta, have been carried out.
Based on the results of the study, several conclusions can be formulated as follows.
The results of the study in Yogyakarta city show that the Uniform Hazard Spectrum (UHS) as a target spectra due to the Shallow crustal is much larger than the Megathrust earthquake source.This result was obtained because the dominant earthquake was a result of the Opak river fault activity located only ± 10 km from the city of Yogyakarta.The Megathrust earthquakes even though they have a larger magnitude but are very far away (> 300 km) of the town of Yogyakarta so that its effect on the hazards is already relatively small.The directivity factor D f and Risk targeted Factor R f respectively to calculate the maximum spectral demand and to ensure the occurrence of uniform risk of building collapse has been used.The value of directivity factor D f > 1.0 [30] , while according to The Indonesian Seismic Loading Code, the risk targeted factor in the short and long period at Yogyakarta city has a value of R f < 1.0.With respect to the uniform hazard spectrum, risk targeted spectrum demand is greater by the average of 12.3% for the entire period T.
Results of the deaggregation process indicated that the dominant magnitude and the distance that contributed the most to the PGA at Yogyakarta city was M D = 6.5 and R D = 14.5 km.This value is consistent with the characteristics and location of the Opak river fault as the primary earthquake source, which is very close to the city of Yogyakarta.Earthquakes due to fault activities are generally much smaller than the megathrust earthquake in the subduction zone.Based on M D and R D , 12-SGMs have been obtained where its actual spectra have matched to the target spectra.The 12 SGMs are relatively varied in both PGA and its A/V ratio.Since the spectral target at bedrock level has a high-frequency content, then almost all of the chosen 12 SGMs also have high-frequency content.The results of spectral matching and the testing of the earthquake significant duration D 595 of the 12 earthquakes have met the requirements.The 12 SGMs have PGA 0.183 -0.276g with an average of 0.232g, and the A/V ratio = 0.771 -2.117, with an average of 1.329, has been found.The SGM significant duration t e = D 595 = 13.04-24.63s(with an average of 15.92s) still falls in the range of the value, according to Eq. 4) i.e between 10.53-25.95s.

Fig. 5
Fig. 5 is a percentage contribution to hazard resulting from the deaggregation process computed according to Eqs. 2 and 3.As shown in the figure, the UHS deaggregation is based only on the shallow crustal earthquake source.The results of the deaggregation

Table 1 .
Source model of subduction earthquake mechanism.
Fig. 3. Earthquake source mechanism of the YSP.

Table 3 .
Prospective of synthetic ground motions