The waveform inversion of mainshock and aftershock data of the 2006 M6.3 Yogyakarta earthquake

This study comprehensively investigates the source mechanisms associated with the mainshock and aftershocks of the Mw = 6.3 Yogyakarta earthquake which occurred on May 27, 2006. The process involved using moment tensor inversion to determine the fault plane parameters and joint inversion which were further applied to understand the spatial and temporal slip distributions during the earthquake. Moreover, coseismal slip distribution was overlaid with the relocated aftershock distribution to determine the stress field variations around the tectonic area. Meanwhile, the moment tensor inversion made use of near-field data and its Green’s function was calculated using the extended reflectivity method while the joint inversion used near-field and teleseismic body wave data which were computed using the Kikuchi and Kanamori methods. These data were filtered through a trial-and-error method using a bandpass filter with frequency pairs and velocity models from several previous studies. Furthermore, the Akaike Bayesian Information Criterion (ABIC) method was applied to obtain more stable inversion results and different fault types were discovered. Strike–slip and dip-normal were recorded for the mainshock and similar types were recorded for the 8th aftershock while the 9th and 16th June were strike slips. However, the fault slip distribution from the joint inversion showed two asperities. The maximum slip was 0.78 m with the first asperity observed at 10 km south/north of the mainshock hypocenter. The source parameters discovered include total seismic moment M0 = 0.4311E + 19 (Nm) or Mw = 6.4 with a depth of 12 km and a duration of 28 s. The slip distribution overlaid with the aftershock distribution showed the tendency of the aftershock to occur around the asperities zone while a normal oblique focus mechanism was found using the joint inversion.


Introduction
In the early hours of , at approximately 05.54 local times, an earthquake struck the city of Yogyakarta. The information from the Agency for Meteorology, Climatology, and Geophysics of Indonesia (BMKG) showed it had a magnitude of Mw = 6.3, centered at 8.003E and 110.320 S with a depth of 11.87 km at the southeast of the city as shown in Fig. 1. The event was declared the deadliest shallow earthquake in the country (BMKG 2006) and destroyed most of the infrastructures in the region of Yogyakarta and Klaten, Central Java Province. The reports by the National Development Planning Agency (BAPPENAS) in cooperation with the Yogyakarta Special Region Government (DIY), Central Java Provincial Government, and international partners in 2006 showed the earthquake killed more than 5000 people, with 38.000 injured, 423.000 evacuated, and 156.000 buildings destroyed (Bappenas 2006).
The mechanism used to obtain valuable information on the earthquake source parameters such as magnitude, fault orientation, stress drop, and source process varied with Open Access *Correspondence: ws@ugm.ac.id 1 Geoscience Research Group, Department of Physics, Faculty of Mathematics and Natural Sciences, Universitas Gadjah Mada, 55281 Yogyakarta, Indonesia Full list of author information is available at the end of the article different agencies such as Indonesia Meteorology, Climatology, and Geophysics (BMKG), United States Geological Survey (USGS), National Earthquake Information Center (NEIC), Kandilli Observatory and Earthquake Research Institute (KOERI), European Mediterranean Seismological Center (EMSC), and Institute de Physique du Globe de Paris (IPGP).
Several debates on the source are also in progress with most researchers reported to have believed the earthquake did not originate from the geological fault along the Opak River due to the distribution of the aftershock towards 10-15 km on the east side (Fukuoka et al. 2006;Walter et al. 2007;Wulandari et al. 2018). Some others also believe some minor faults were reactivated by the earthquake on the east side of the Opak fault (Diambama et al. 2019;Budiman et al. 2019;Irham et al. 2014). According to Saputra et al. (2018), the Opak River consists of 56 faults with a maximum displacement of 2.93 while Nakano et al. (2006) found a different focal mechanism as shown in Fig. 1. Moreover, several factors have been indicated by researchers to be influencing this difference, such as the lack of seismometer networks of BMKG which has been existing before 2006, thereby, causing variations in the level of location accuracy, source depth, and fault orientation (Ma and Eaton 2011;Saunders et al. 2016).
The parameters discovered by several institutions to indicate the source of the earthquake are presented in Table 1. The majority were observed to have concluded the epicenter is on the east side of the Opak fault and this also indicates consistency based on the aftershock distribution (Walter et al. 2008). Kawazoe and Koketsu (2010) determined the focal mechanism using the waveform inversion method and compared the findings with the data observed using the CRUST 2.0 velocity model (Bassin et al. 2000) as well as the left-lateral strike-slip and reverses dip-slip. Meanwhile, Nakano et al. (2006) examined the focal mechanism dominated by strike-slip and horizontal coseismal displacement was discovered to have occurred around Bantul and Yogyakarta in the south and southwest direction with the focal mechanism found to be in the form of left-lateral strike. Similar results were also obtained by Fig. 1 Focal mechanisms results from several agencies. The yellow star represents the location of the epicenter, the black and white shoreline indicates the focus mechanism, the dark blue line is the Opak fault, and the light blue line is the minor fault scattered east of the Opak fault Saputra et al. Geosci. Lett. (2021)  This present study was, therefore, conducted for two reasons. The first involved using moment tensor inversion to determine the type of expansion experienced based on the mainshock and three aftershocks Mw > 4. The second was to determine the earthquake source model by understanding the distribution of cosmic slips in space and time using the joint inversion method. It is important to note that the moment tensor and joint inversions were calculated based on the fit between the observed and synthetic waveforms with the smallest variance.

Moment tensor inversion mainshock and aftershocks Data
The data used for the mainshock moment tensor inversion were obtained from the Incorporated Research Agency for Seismology (IRIS) while those for the aftershocks were retrieved from the Geofon and European Integrated Data Archive (EIDA) with the focus on the 8th, 9th, and 16th of June 2006 through (http:// ds. iris. edu) and (http:// eida. gfz-potsd am. de/ webdc3), respectively. The distribution of the IRIS and BMKG networks for all stations is, however, presented in Fig. 2. The waveform data used to distribute the aftershocks were from the non-permanent seismometer installed one week after the mainshock around the Opak fault from 3 to 7 June 2006 by the Universitas Gadjah Mada in cooperation with the GFZ, Potsdam, Germany (Wulandari et al. 2018). It consists of ten stations with three-component sensors and the experiment found approximately 524 aftershocks events which were relocated and plotted to visualize the spatial distribution as shown in Fig. 2. Meanwhile, the aftershocks data plotted were retrieved from Wulandari et al. (2018). It is important to note that the mainshock data were relocated using the SeisComP3 program while the aftershocks distribution data were through the Non-LinLoc method (Wulandari et al. 2018).

Method
The waveform was corrected from its instrument response using the Fortran code developed by Yagi (2006) and the mainshock waveform was windowed with a data length of 100 s which started 5 s before the P-wave arrival time. Meanwhile, the aftershock waveform was windowed with a length of 50 s starting 5 s before the P-wave arrival time. The filters used during the process were in accordance with those applied in previous research as presented in Tables 2 and 3. Moreover, the filter and velocity model were both selected based on the smallest variance value. The variance is, however, the comparison between the observed and calculated synthetic waveform using Eq. (1) (Yamanaka and Ishida 1996;Ito et al. 2004).
where u obs j is the observed waveform at each j station and u cal j is the calculated waveform at each j station at time t. The optimal inversion model was obtained by computing the moment tensor inversion, which iterates at appropriate resampling time values from 0.1 to 10 s.
A grid search method ranging from 1 to 25 km was applied in this study to determine the depth of the hypocenter for both the mainshock and aftershock in each frequency range with a bandpass filter. The waveform was converted from velocity to displacement using a sampling time of 0.2 s-10 s and the Green function calculated using five velocity models which are AK135 (Kennett et al. 1995), CRUST2.0 (Bassin et al. 2000), PREM (Dziewonski and Anderson 1981), JB (Jeffreys and Bullen 1940), and Koulakov (Koulakov et al. ( where G jq is the Green's function of time t and τ is a unit step from the source positions of x, y, and z . M q is the moment tensor element, V is the volume of the (2)   earthquake source space, q is the number of free components for the second pair selected, and e 0 is the observation error. Meanwhile, the focal mechanism represents a point source model. Therefore, Eq.
(2) can be simplified into a vector shape to form Eq. (3) as follows: It is possible to rewrite this in vector form as follows: where T (t) is the source time function in the source centroid, x c ,y c ,z c , d, and e j are data vector errors with N-dimensions, m represents the five dimensions of the model vector parameter, and G is the N × 5 coefficient matrix. In Eq. 4, u ud and u ns are the observation waveform data for vertical and horizontal components while (t 1 ), (t 2 ) . . . .(t n ) are the arrival times of the observed waveforms at a station.
The solution to the matrix equation was determined using the least square method when the observed waveform and the convolution of the Green's function with the time function of the source are known. Meanwhile, Green functions for the near-field data were calculated using the extended reflectivity method developed by Kohketsu (1995), Kuge (2003) and Yagi (2006). The source time function was also used to determine the moment tensor solution while the slip-rate function was applied to the waveform data with a simple triangular function.

Joint inversion Data
The wave data for P and S were downloaded in the form of 16 teleseismic waveform body with the vertical components recorded on the IRIS-DMC station which was downloaded via http:// ds. iris. edu. Moreover, the teleseismic data for the mainshock event on 27 May 2006 were downloaded from The International Federation of Digital Seismograph Networks (FDSN) on January 2, 2019. The selected data have a high signal-to-noise ratio and covers a good azimuth range with the source and station distances used observed to be within the range of 30°-90°. The teleseismic station configuration is, however, indicated in Fig. 3a. Meanwhile, the near-field data which  Fig. 3b. The teleseismic data used for the wave inversion calculation consist of 12 stations, windowed with a data length of 90 s, 5 s before the P wave, filtered with a bandpass filter in the frequency range of 0.008-0.1 Hz, and converted to displacement with a sampling time of 0.5 s. Meanwhile, the near-field data used 14 components from 5 broadband seismometer stations including four which are TNG, BJI, KMMI, and DNP from the BMKG network and one which is XMIS from the IRIS-DMC network. The data were windowed with a data length of 180 s, 5 s before the P wave, filtered using a bandpass filter with a frequency range of 0.05-0.15 Hz, and converted to displacement with a sampling time of 0.5 s. The first P wave was, however, manually picked with the help of the Seismic Analysis Code (SAC) program when it arrived at each seismogram (Helffrich et al. 2013).

Method
The rupture process of the Yogyakarta earthquake source was determined using the equation from Yagi et al. (2004) which was developed by Yoshida and Koketsu (1992). The equation represents the rupture process as a temporal and spatial slip distribution in the fault plane and the wave received by the station is the observed waveform at each station j given as indicated in Eq. 5: where W obs j is the observed waveform at station j , X mnkl is the slip component at each mn sub faults and each l th time change, g mnkj (t) is the basic seismic wave Green function of the point source with slip units on each mn sub faults, τ is the instantaneous time at the subfault points, T mn is the basic function of the start time for each sub faults, and e j is the gaussian error with the variance σ j . Meanwhile, the Green function used for teleseismic data was calculated using the Kikuchi and Kanamori (1991) method while the seismic velocity model applied to calculate Green data near-field and teleseismic functions was Koulakov and Jeffreys Bullen (JB), respectively, as indicated in Table 3.
A numerical method was applied to the standard waveform inversion scheme to obtain a more objective source model (Harzell and Heaton 1983;Yoshida 1992). Moreover, the inversion code developed by Yagi et al. (2004) was also adopted to calculate inversion due to its ability to represent rupture as a slip distribution over time and space in a fault plane. The initial model of the fault plane was determined from the area distribution while Papazachos et al. (2004) was also applied to calculate the slip lengths and areas. Furthermore, the fault plane was divided into MxN sub faults with the length and width represented by d x and d y , respectively. The slip-rate function on each sub-default was determined using a simple triangular series of functions for the length of the dislocation or rise time, τ . Therefore, Eq. (5) was simplified into a vector form as follows: Fig. 3 Near-field data station configuration and body wave teleseismic data. a The station configuration for the body wave teleseismic data is provided as a red inverted triangle symbol using 12 stations and the distance of the teleseismic data station to the source is 300-900. The yellow star in the pink box is the epicenter position and the near-field station configuration. b The near-field data station is marked with a red inverted triangle symbol and the yellow star is the epicenter of the Yogyakarta earthquake mainshock while the station-to-source distance for the near-field data is 0°-5. 5° Saputra et al. Geosci. Lett. (2021) 8:9 where A j is a matrix having dimensions in the form of some data points at each station j× the number of model parameters N j .
Changing the number of model parameters leads to an unstable solution due to the ability of a small change in data to cause a big difference. Therefore, smooth limits were applied to the slip distribution over time and space to produce a stable inversion.
There are two constraints, spatial and temporal, discovered in the seismic wave analysis conducted using observational data. This method was, however, used to obtain the best model parameters calculation by minimizing the number of residual squared, S, which was determined using the following equation: where T is the matrix of N 1 xN a (N 1 = MNLK ) , D is the matrix of N 2 xN a (N 2 = MNK ) , MNK is slip component K th on the sub faults MN th and timestep L th , while N a is the model parameter.
The relationship of Eq. 5 coded by Yagi et al. (2004) showed the smoothing constraints of space and time as provided in the following equation: where a klmn is the model parameter determined from the observed data, X kmn (ξ ) is the basic slip function of slip k th , and T lmn (t) is the time base function of the l th slip time. Meanwhile, the smoothing constraint with respect to time was obtained by simplifying Eq. 8 into Eq. 9: where X mnk0 = X mnkl = 0 . This was further converted into a vector form: where T is the matrix of N 1 × N a (N 1 = MNLK ) . Meanwhile, the distribution of slip to space for the smoothing constraint is presented in the following vector form: where D is the matrix of N 2 × N a (N 2 = MNK ).
The least-square method was applied to determine the best estimate for the model parameters by providing the covariance value σ (Cheng et al. 1985). Meanwhile, it is impossible to directly determine the values of σ t and σ d ∂ 2 ∂t 2 a klmn X kmn (ξ )T lmn (t) + e j = 0, (9) X mnkl(l−1) − 2X mnkl + X mnk(l−1) + e j = 0, (10) Tx + e j = 0, (11) Dx + e j = 0, while σ j can be selected from the data quality. Therefore, σ t and σ d were selected objectively using Akaike's Bayesian Information Criterion (ABIC) method (Akaike 1980) to obtain the optimal values with the grid search method. The same equation was also used in Fukahata et al. (2003Fukahata et al. ( , 2004 and it is usually represented as follows: where N is the total number of observed equations and C is a constant. The grid search method was applied in this study to obtain the optimal values of σ t and σ d while the non-negative least-squares (NNLS) method was used to solve the least-squares problem with a positive constraint on the model parameters.

Moment tensor inversion of mainshock event
The mainshock source mechanism was determined by the grid search method using different frequency ranges and velocity models as shown in Tables 2 and 4 and the results are presented in Fig. 4. Meanwhile, the velocity waveform was converted to displacement using a sampling time of 0.1-10 s in order to select the appropriate stable result. Figure 4 shows the beachball variations and the time sampling with the smallest variance value for each velocity model was found to be 0.5 s while the frequency was 0.01 -0.05 Hz. This sampling time was observed to be very stable with the five velocity models while those above 1 were precarious. Moreover, the frequencies with low variance were found at low values between 0.01 to 0.05 Hz. This means the appropriate velocity model to describe the subsurface structure of the study area is Koulakov with a sampling time of 0.5 s and a frequency of 0.01 to 0.05 Hz. The minimum variance value for the grid search method was 0.234 based on the comparison between the observed and synthetic waveforms as shown in Fig. 5. Furthermore, the parameters of the Yogyakarta seismic source on May 27, 2006, include the moment 0.2808 E + 19 (Nm) which is equivalent to the moment magnitude M w = 6.2. Meanwhile, the fault parameters obtained for Nodal 1 were strike 353.7°, dip 43.4°, and slip − 126° while Nodal 2 had 218.7°, 56.2°, and − 61°, respectively, and the hypocenter was 12 km.
The results of the moment tensor inversion from nearfield data for the mainshock are shown by Nodal 2 with the type of faulting found to be strike-slip and dipnormal. The selection of Nodal 2 was, therefore, based on the Opak fault lithology and geological maps with a (12) Saputra et al. Geosci. Lett. (2021) 8:9 southern direction. This, therefore, made strike 218.7° the most suitable strike direction. The faulting type is, however, depicted as a beachball as shown in Fig. 6.

Moment tensor inversion of aftershocks event
The aftershock mechanism was determined in the same way as the mainshock by using the grid search method as observed in previous studies at different velocity models and frequency ranges as indicated in Tables 3 and 4, respectively. The aftershock reversal results for June 8, 9, and 16 are presented in Figs. 7, 8, and 9, respectively. Meanwhile, the waveform velocity was converted to displacement using a sampling time of 0.2 to 10 s.
The graphs in the figures are presented as a function of time sampling and variance for each model while the speed and range of the bandpass filter are indicated in Table 2. The sampling time of 1 s and above for the aftershock on June 8, 9, and 16 in the Fig. 4 The grid search method to obtain the best solution for the mainshock focal mechanism as a function of time sampling and variance for each speed and range for the bandpass filter models as shown in Table 2. The blue, yellow, red, brown, and green color beachballs indicate the filtering used by Nakano, Suardi, Mikumo, Yagi, and this study, respectively. The velocity models used include a AK135, b CRUST, c PREM, d JB and e Koulakov Saputra et al. Geosci. Lett. (2021) 8:9 five velocity models have the same trend which is the focal ball and unstable variance values. Meanwhile, the 0.5-s sampling time for all the aftershocks was very stable, therefore, the Koulakov velocity model was applied because it is very representative of the subsurface structure of the study area. The lowest dispersion frequency recorded for the five velocity models and aftershocks was 0.1-0.3 Hz while the smallest variance values for 8, 9, and 16 were 0.2178, 0.1949, and 0.2432, respectively. The parameters of the inversion result for the aftershock source are presented in Table 5. Figure 10 shows the fitting shape of the observed waveform and the synthetic waveform displacement for the aftershocks and the waveform fitting produced is the best result from the grid search process.
The variation of the fault plane types for the aftershocks is presented in Table 5 and the typical dip-normal was found to have been used on June 8 while strike-slip dominated June 9 and 16. Moreover, the magnitude has a similar result with the catalog from International Seismological Center (ISC) which was 4.4, 4.1, and 4.0 for 8, 9, and 16, respectively, as shown in Fig. 11. Meanwhile, the variance value obtained from the inversion is also presented in Table 5.  5 The fitting of the displacement between the observed and synthetic waveform from the Yogyakarta earthquake mainshock recorded by BJI, JCJ, TNG, SBJI, KMMI, and XMIS stations. The blue and red lines are the observed and synthetic waveforms, respectively, and the value above the waveform is the maximum amplitude of the observed waveform Saputra et al. Geosci. Lett. (2021) 8:9 The aftershock moment tensor inversion shown in Table 5 serves as the basis to examine the consistency of the fault plane caused by the mainshock. The beachball solution from this inversion is, however, indicated in Fig. 11. Figure 11 shows the inversion result for the mainshock and aftershocks moment tensor and the mainshock focal mechanism was found to be a dip-normal strikeslip type of fault and the same was observed for June 8 while June 9 and 16 had a strike-slip. There is, therefore, a strong suspicion that June 9 and 16 aftershocks originated from a different fault plane than the mainshock with the report presented by USGS showing strike-slip was predominant.

Joint inversion
Akaike's Bayesian Information Criterion (ABIC) was adopted to obtain more objective parameters (Akaike 1980). The Green function for teleseismic body waves was calculated using the Kikuchi and Kanamori (1991) method while the technique developed by Kohketsu (1985) was used for the near-field data.

Fault model
The fault model used in this study was the Hazkel finite model, which assumes there is no change in the slip angle along the rupture and that faulting occurs in a single fault plane. This means fault dimensions are essential to determine the detailed global features of the fault plane and coseismic rupture of an earthquake. The seismic traces in the fault plane can, however, be used to estimate its extent. Therefore, mapping the distribution of the Yogyakarta earthquake aftershock from several previous studies shows the area of the Opak fault plane is 28 km long and 15 km wide (Walter et al. 2008;Anggraini 2013). Meanwhile, Papazachos et al. 's (2004) equation showed the length is 30 km and the width is 20 km. The initial model of the fracture area used in this study was divided into stage 1 and stage 2 to obtain the area of the fault plane which corresponds to the actual area of the Opak Fault. Stage 1 was based on the aftershocks distribution while Stage 2 was based on the equations of Wells and Coppersmith (Wells and Coppersmith 1999). The fault plane on stage 1 with a strike length of 32.5 km and a dip width of 15 km was divided into 15 × 6 sub faults Fig. 6 Focal mechanism of the inversion model using Koulakov's velocity model, which is illustrated by a green color beach ball Saputra et al. Geosci. Lett. (2021) 8:9 with each having a length of 2.5 km. The dimensions of the fault plane on stage 2 are smaller than stage 1 with a long strike of 17.5 km and a dip width of 15 km, which was further divided into 7 × 6 sub-faults. Moreover, the hypocenter position used was the result of the relocation from the mainshock released by USGS using the software sesicomp3. The fault plane model is presented in Fig. 12. Moreover, the velocity model used to calculate Fig. 7 The grid search method to obtain the best solution for the focal aftershock mechanism on June 8 as a function of time sampling and variance for each model speed and bandpass filter range as shown in Table 2. The blue, yellow, red, brown, and green color beachballs indicate the filtering used by Pinar, Kasemsak, Miyatake, and this study, respectively. The velocity models used include a AK135, b CRUST, c PREM, d JB and e Koulakov Saputra et al. Geosci. Lett. (2021) 8:9 the teleseismic body wave data was Jeffreys Bullen (JB) while five velocity models were used to calculate the near-field data as shown in Table 4. Meanwhile, the JB velocity model for the teleseismic data was changed into four near-field data rate models to obtain the smallest variation. Fig. 8 The grid search method to obtain the best solution for the focal mechanism of the June 9 aftershock as a function of time sampling and variance for each of the velocity models and bandpass filter ranges as shown in Table 2. The blue, yellow, red, brown, and green color beachballs indicate the filtering used by Nakano, Suardi, Mikumo, Yagi, and this study, respectively. The velocity models used include a AK135, b CRUST, c PREM, d JB and e Koulakov Saputra et al. Geosci. Lett. (2021) 8:9 The rupture velocity was set to be equal with the maximum shear wave velocity of 4.6 km/s used in Koulakov et al. in order to save computation time. The waveform data applied to implement the slip-rate function was a simple triangular function sequence with a rise time of 1 s while the dislocation time and rupture propagation Fig. 9 The grid search method to obtain the best solution for the focal mechanism of the June 16 aftershock as a function of time sampling and variance for each of the velocity models and bandpass filter ranges as shown in Table 2. The blue, yellow, red, brown, and green color beachballs indicate the filtering used by Nakano, Suardi, Mikumo, Yagi, and this study, respectively. The velocity models used include a AK135, b CRUST, c PREM, d JB and e Koulakov for each sub faults were used to determine the slip-rate function duration. Moreover, the hypocenter depth was determined by calculating the mainshock moment tensor inversion using the grid search method. The fault plane parameter and hypocenter depth were further used as the initial input parameters to calculate the joint inversion.

Source rupture process
The rupture of the Yogyakarta earthquake source on May 27, 2006, was due to a joint inversion as indicated in Fig. 13, which shows the source mechanism, moment rate function, and the coseismic slip distribution. The focal mechanism from the joint inversion is similar to the initial input with the fault parameters recorded to be strike 217.6°, dip 56°, and slip − 53° as shown in Fig. 13a while the source duration of moment rate function length was 28.5 s as indicated in Fig. 13b. Bilek and Lay (2002) have previously stated that shallow earthquakes close to subduction zones have weak geological material because they are composed of young sedimentary rocks (Pribadi et al. 2014). Moreover, Saputra et al. (2018) also reported geological conditions in the eastern and western parts of the Opak fault to have different structures. The east has a more slightly fragile structure, which makes it susceptible to landslides while the west has a young Merapi sedimentary formation structure and this difference leads to a variation in the P and S wave velocity, thereby, making the rupture duration of the source reach 28.5 s. It is important to note that the total seismic moment produced was M0 = 0.4311E + 19 (Nm) or Mw = 6.4 as shown in Fig. 13b and this is similar to the M0 = 0.422E + 25 (Nm) or M0 = 6.3 released by USGS. Meanwhile, the energy moment was released at the 6th second after the first initial break and the total slip vector produced was 0.78 m as shown in Fig. 13d. The slip moved in an NS direction from the hypocenter point, later to the SW direction towards the surface, and moved in a dip direction at the 15th second to the end. This shows the Yogyakarta earthquake fault type is more dominant dip-normal than strike-slip as shown in Fig. 13c. Figure 14a shows a series of snapshots showing the Yogyakarta earthquake dislocation distribution to the surface area from 0 to 28 s. The change in dislocation was observed not to have occurred at the initial break but at the third second by moving southward as far as 10 km up to the eighth second after which it moved northward for 20 km up to the twelfth second and the slip movement was seen moving in the direction of 15 km from the fifteenth up to the last second. The maximum slip recorded was 0.72 m at the 10 km north of the hypocenter. This means the first and second asperity zones are both 10 km south of the hypocenter as shown in Fig. 13c. Meanwhile, the slip rate of snapshots per second is illustrated in Fig. 14b and the average rupture speed was found to be 4 m/s. Figure 15 compares the observed waveform in blue with the synthetic waveform in red displacement for the near-field data in (a) and body wave teleseismic data in (b). The variance value was recorded to be 0.3412 and this is quite large, even from the waveform fitting perspective, and observed to be quite good at the beginning of the P wave arrival. The limited quality of the signal recorded by the teleseismic station is also a factor in the variance value. This is due to the fact that earthquakes with 6.3 magnitudes are less well recorded at long-range stations and also affected by layer inhomogeneity. Moreover, limited azimuthal coverage of near-field data also has effects on the quality of inversion calculation. It is important to note that the near-field and teleseismic data do not have very good quality, but showed a fairly good resolution quality and this is one of the advantages of simultaneous inversion using ABIC optimization (Akaike 1980). Therefore, the variance values obtained are the best results from the joint inversion.   Figure 16 illustrates the Yogyakarta earthquake coseismic slip overlaid with the distribution of the relocated aftershock. It also explains the aftershock distribution pattern with the coseismic slip distribution. Figure 17 shows the depth profile of the aftershock is divided into a few slices which are A-Aʹ, B-Bʹ, and C-Cʹ, and the most dominant aftershocks were observed to have occurred at a depth of 10 km as indicated in Fig. 17b. Meanwhile, the C-Cʹ section shows an empty space with no energy released as shown in the slip distribution of Fig. 16 and the aftershock was discovered to have lain on a high enough slip. This pattern shows the faults are in the same plane, but those located on the west are different. Moreover, the aftershock moment tensor inversion previously produced on June 9 also indicated different types of faults and this further shows the complexity of the faults in the study area.

Discussion
The moment tensor inversion obtained using near-field data showed the mainshock fault source parameters of the Yogyakarta earthquake which occurred on May 27, 2006 were strike 218.7°, dip 56.2°, and slip − 61° and discovered to be caused by an oblique fault. This is in line with the findings of previous studies that the Yogyakarta mainshock was not caused by pure strike-slip but has a dip-normal component. Moreover, Kawazoe and Koketsu (2010) divided earthquake source mechanism into two types which are the left-lateral strike-slip faulting and reverse dip-slip faulting and this further indicates the possibilities of two different fault segments for the source mechanism. Tsuji et al. (2009), however, found a movement in the direction of the fault line on the east side along the Opak fault line. These results confirmed the joint inversion findings indicated in the slip distribution of Fig. 11c which showed a slip shift towards the dip at the 15th second. The most significant slip also occurred at the second asperity at 0.72 m.
The moment tensor inversion for June 8, 9, and 16 aftershocks were calculated using temporary near-field data and the fault source parameters also showed the complexity of the Opak fault (Budiman et al. 2019;Saputra et al. 2018;Anggraini 2013). The aftershock of June 8 was observed to have dip-normal faults while 9 and 16 had strike-slip as shown in Table 5. This means the Yogyakarta earthquake mainshock activated several active faults on the eastern side of the Opak fault (Anggraini 2013;Budiman et al. 2019;Irham et al. 2014).
The rupture process of the earthquake source was also determined using joint inversion and the earthquake source mechanism, moment rate function, and the coseismic slip distribution were discovered. The focal mechanism from the inversion was observed to be similar to mainshock moment tensor inversion with the fault parameters recorded to be strike 217.6°, dip 56°, and slip − 53°. This, therefore, shows the consistency between the moment tensor and joint inversions. Moreover, a

Fig. 11
The mainshock and aftershock focal mechanism for the Yogyakarta earthquake on May 27 2006 is represented in green while the red color small circle indicates the distribution Fig. 12 The assumed fault plane used in the joint inversion calculation. The black line is the assumption at stage 1 while the green line is the assumption at stage 2. Yellow stars and red circles are the distribution of epicenter locations and aftershock relocations while the green beachball is the focal mechanism of the inversion mainshock dip-normal fault type was recorded in the joint inversion with a smaller slip than the moment tensor inversion. Yagi (2006) in Ohsumi and Baba (2007) showed the rupture due to the Yogyakarta earthquake was in two stages with the first stage observed to have the first subevent in the southwest to northeast directions while the second stage had the second subevent. The findings of this study also indicated two weak zones which were located Fig. 13 The joint inversion of the Yogyakarta earthquake on May 27, 2006. a The focal mechanism of the Yogyakarta earthquake. b Functions of the moment rate and source duration. c The distribution of the coseismic slip vector with the yellow star used to represent the hypocenter position, the color of the contour represents the size of the slip dislocation, and the white dots are the position of the Asperity zone. d The graph of the dislocation in the asperity zone against time Fig. 14 Series of snapshots projected orthogonally onto the surface plane. a Distribution of dislocations per second from 0 to 28 s. b Distribution of the slip rate per second from 0 to 28 s. The yellow star is the epicenter of the mainshock 10 km south of the hypocenter as the first Asperity and 10 km to the north as the second Asperity. Kawazoe and Koketsu (2010) also identified two zones of violence near the hypocenter and in the southwest and this raised the suspicion that the rupture was from two different events. These results are also confirmed by the aftershock moment tensor inversion with the type fault recorded for June 9 and 16 aftershocks different from June 8. This means the reactivation of several minor faults possibly occurred on the east side of the Opak fault (Anggraini 2013;Walter et al. 2008;Saputra et al. 2018;Budiman et al. 2018).
The coseismic slip distribution overlaid with the aftershock relocation distribution as indicated in Fig. 15 showed the aftershock spread out parallel to the Opak fault. Meanwhile, another pattern was discovered to the west of the Opak fault distribution and believed to have originated from different fault segments. This, therefore, confirms the findings from the moment tensor inversion for June 9 and 16 aftershocks. Fig. 15 Comparison of the observed and synthetic waveforms. The blue is the observed waveform while the red is the synthetic waveform. The value above the waveform is the maximum amplitude value of each waveform. a Near-field data. b Teleseismic body wave Saputra et al. Geosci. Lett. (2021) 8:9

Conclusions
The moment tensor inversion of the mainshocks and aftershocks as well as the joint inversions were calculated to obtain a more comprehensive source model for the Yogyakarta May 27, 2006 earthquake. The exact source mechanism was determined using the mainshock waveform inversion after which the source parameters which include dip-normal and strike-slip were used as the input to calculate the joint inversion. Moreover, the ABIC method was adopted to obtain a stable result in the joint inversion (Akaike 1980) and determine the coseismic slip distribution in the rupture area which was found to be bilateral. The slip distribution direction was generally discovered to have propagated along the fault plane in the strike direction and later moved and stopped at the dip direction. It spread from the hypocenter point to the south for 10 km and later moved north for 20 km after which it moved towards the dip-normal. The earthquake was found to have originated from the Opak fault and also reactivated a minor fault on the east side. It was also strongly