1 Introduction

Merapi volcano is located in central Java and surrounded by densely populated areas so that its eruption effect will be hazardous for the people living around the volcano. Merapi has a relatively very short return period of eruptions, i.e., from 4 to 6 years with volcanic explosivity index (VEI) of 1–3 (Andreastuti et al. 2000; Voight et al. 2000). The eruption type of Merapi is dominated by pyroclastic flow caused by dome collapse. High VEIs (more than 4) with an explosive type occur in a longer return period, for instance the 2010 eruption event. Such an explosive type of eruption occurs every ±100 years (Surono et al. 2012). This eruption type has been controlled by deep magma sources with different mechanism (Costa et al. 2013). Based on the result of a petrological study (Costa et al. 2013), there are two magma bodies beneath the Merapi volcano, i.e., at ~30 and ~13 km depths. The DOMERAPI project was conducted partly to understand the internal structure and characteristics of the Merapi volcano.

Many earth scientists have conducted research on Merapi using various methods to understand the process related to physical, chemical and geological parameters beneath the volcano and its summit. Seismological approach can also be used to understand physical parameters beneath a volcano. For example, studies using MERAMEX (Merapi Amphibious Experiment) data have resulted in some publications explaining the relation between subduction zone and volcanic arc in central Java (Wagner et al. 2007; Koulakov et al. 2007; Rohadi et al. 2013; Haberland et al. 2014). These studies depict clearly a low velocity ascending from the subducted slab to the volcanic arc, but they do not show a detailed velocity structure of Merapi due to the lack of data coverage beneath the volcano.

Detailed studies of Merapi have been performed using volcano-tectonic (VT) earthquake data to analyze seismicity and focal mechanism of events beneath the volcano (Ratdomopurbo and Poupinet 2000; Hidayati et al. 2008; Budi-Santoso et al. 2013). They show that the maximum depth of VT earthquakes is about 5 km depth from the summit, so structures below this depth cannot be delineated.

In this study, the data from the BMKG network in the same period of the DOMERAPI project were incorporated in the hypocenter determination and relocation to maximize azimuthal coverage. Epicenters from the DOMERAPI network have an azimuthal gap average of more than 300°. The joint networks have successfully minimized the azimuthal gap to <250° (Ramdhan et al. 2015).

2 Data and method

The DOMERAPI network covers Merapi and Merbabu volcanoes and geological features, such as the southern mountain range of central Java, part of the Kendeng basin and Opak fault (Fig. 1). The network consists of 53 sites of broad-band seismometers operated from October 2013 to mid-April 2015. The network consists of 46 seismometers at the beginning, and then seven seismometers are relocated to different sites to obtain better S/N ratio. The seismometer spacing average is about 4 km. This network design is potential to enhance seismic image resolution in the study area in unprecedented detail. A preliminary work on seismic ambient noise tomography was applied using the network to delineate subsurface structure beneath Merapi and Merbabu volcanoes (Trichandi et al. 2015).

Fig. 1
figure 1

Distribution of seismographic stations of the DOMERAPI network. The inset picture shows the boundary of DOMERAPI network area. The geological map was modified from Smyth et al. (2008). The 2006 earthquake location was taken from global CMT catalogue (Ekström et al. 2016)

Full size image

Broad-band seismometers detect wide frequency signals, so they detect not only local events but also regional and teleseismic events as well. For this study, we selected regional events in the radius of 500 km from the Merapi volcano. Figure 2a, b shows examples of 3-component waveforms from local and regional events recorded by the DOMERAPI network. We identified events and picked P- and S-wave arrival times carefully. Events were determined by at least six DOMERAPI stations for P and S phases or P phase only. The P phase has much better onset quality than S leading to a less number of S phase than P.

Fig. 2
figure 2figure 2figure 2

a Examples of 3-component waveforms from local events recorded by the DOMERAPI network. b The same as a, but from regional event. c Shows the epicentral distribution of earthquakes and stations recording these two events. The left and right figures are close up the study region to show the station distribution recording the local and regional events, respectively

Full size image

The earthquake hypocenters were first determined using the Geiger’s method (Geiger 1912) implemented with the hypoellipse software (Lahr 1999). The method was also applied to determine earthquakes along Cimandiri and Baribis faults in West Java based on BMKG data (Supendi and Nugraha 2016). We relocated the hypocenters using a double-difference earthquake location algorithm (Waldhauser and Ellsworth 2000; Waldhauser 2001). This method was successfully applied to relocate hypocenters in the Indonesia region based on data from BMKG network (Ramdhan and Nugraha 2013; Cahyaningrum et al. 2015; Sabtaji and Nugraha 2015; Utama et al. 2015). The hypocenter relocation was conducted carefully to obtain precise hypocenters before conducting seismic travel time tomography for our future work.

3 Results and discussion

In this study, we incorporated events recorded by the BMKG seismometers outside the DOMERAPI network as shown in Fig. 3 to maximize azimuthal coverage. There are 33 broad-band seismometers selected from the BMKG network. XMIS is one of the Australian seismographic stations that was used in the processing in order to better constrain event locations, especially to the south of Java. The hypocenter determination and relocation hypocenter used a 1-D P velocity model for central Java (Koulakov et al. 2007) and v P/v S ratio of 1.73 obtained from the Wadati diagram constructed using the combined DOMERAPI and BMKG data sets (see Wadati 1933).

Fig. 3
figure 3

a Epicentral distribution from the DOMERAPI and BMKG networks before hypocenter relocation. b The same as a, but after relocation. Symbol yellow triangles depict the location of seismographic stations of the BMKG network used to determine earthquake hypocenters along with the DOMERAPI stations

Full size image

There are 464 earthquakes determined from the data recorded by the DOMERAPI and BMKG networks. After hypocenter relocation, we end up with 399 events. There are some earthquakes that could not be relocated which may be due to poorly clustered events. A comparison between epicenter distributions from the DOMERAPI and BMKG networks before and after relocation is shown in Fig. 3. Most of the resulting hypocenters are located outside the DOMERAPI network, so they can be useful for imaging the deep structure of Merapi in our future work on seismic tomography. Our new data set indicates that there are still many events to the east of the Opak fault events after nearly a decade of the Yogyakarta earthquake in 2006. In addition, a thrust fault to the south of central Java, which may be related to a backthrust, is clearly delineated by the relocated events.

Epicenter shifts after relocation and a rose diagram depicting the number of relocations for 5° azimuth bins are given in Fig. 4. It is shown clearly that most of the shifts are perpendicular to the trench. The average epicenter shift is ~5 km, and the maximum shift is ~27 km. The majority of the epicenter shifts are northward which may be caused by the seismographic station distribution, in which almost all of the stations are located in Java Island. The shifting perpendicular to the trench has been interpreted as most likely due to an effect of a 3-D velocity model incorporating the slab, because teleseismic rays spend a lot of time in the slab (see, e.g., Pesicek et al. 2010). In this study, however, we used a lot of local data and 1-D velocity model, but we still obtained a similar result. So the shifting perpendicular to the trench seems to be caused by the station distribution, which is generally located on land. The horizontal errors (X and Y error) and vertical errors (Z error) decrease significantly after hypocenter relocation. The majority of errors for all direction are <0.3 km. Figure 5 shows the error distribution before and after hypocenter relocation. The result of the hypocenter relocation is used to analyze the seismotectonics in the study region as depicted by the following figures.

Fig. 4
figure 4

a Epicenter shifts owing to event relocations. Rectangles depict the location of vertical cross sections A–A′ and B–B′ shown in Figs. 6 and 8, respectively. b Rose diagrams showing number of relocations for 5° azimuth bins. The shifting values are normalized first to be plotted in the diagram. The maximum value in the rose diagram is 1.0, i.e., a normalized value. The increment scale of 0.2 in the radial direction is shown

Full size image
Fig. 5
figure 5

a Plot of errors in X (east-west) direction before hypocenter relocation (red square) and after hypocenter relocation (blue square). b The same as a, but in Y (north-south) direction. c The same as a, but in Z direction

Full size image

There are two tectonic features detected clearly by the relocated hypocenters. They are: (1) a nearly vertical cluster of events that may be related to a backthrust to the south of the study region and (2) the active Opak fault to the south of Merapi, which is indicated by small events even after a decade of the Yogyakarta earthquake of May 26, 2006, at 22:54 UTC (Ekström et al. 2016).

A vertical section (A-A′) crossing the Opak fault is shown in Fig. 6. The cluster of events beneath the top of slab, from a depth of ~40 km down to ~115 km, may be related to a double seismic zone as discussed by Koulakov et al. (2007), although we could not delineate it clearly. Notice the existence of fixed depths before hypocenter relocation, which was successfully overcome after the hypocenter relocation.

Fig. 6
figure 6

a Vertical cross sections A-A′ showing the distribution of hypocenters a before and b after relocation. Blue line depicts the slab 1.0 model from USGS (Hayes et al. 2012). The location of cross sections is shown in Fig. 4

Full size image

Figure 7 depicts the seismic activity occurring to the east of the Opak fault during the DOMERAPI project period. This fault has been interpreted to be related to the Yogyakarta earthquake occurring on May 26, 2006. The earthquake mechanism is sinistral strike slip with M W 6.3 (Nakano et al. 2006; Ekström et al. 2016). The earthquake occurred to the east of the Opak fault line. The earthquake source may come from a hidden fault because no geomorphic feature was well noticed before the earthquake (Tsuji et al. 2009). Geological condition shows that surrounding area beneath the hidden fault is dominated by sediment covers. Usually an earthquake related to an active fault occurs at the basement part of continental crust. For the case of a thick sediment cover, the slip that is created by an earthquake is not much stronger to develop surface rupture and make a geomorphological feature like surface faults. It is consistent with the geological condition, in which the hypocenter of the earthquake was located beneath an alluvial deposit (Smyth et al. 2008). Based on the relocation result, the depth of the majority of events to the east of the Opak fault is <25 km. This is in agreement with the crustal thickness in the area, i.e., around 35 km (Wölbern and Rümpker 2016). From the hypocenter depths, it is clear that these events indeed occurred in the basement part of continental crust. At this stage, we speculate that the cluster of events parallel to Opak fault is alternatively related to the fault. Our preferred interpretation is that the Opak fault is dipping eastward. However, a further investigation is still needed.

Fig. 7
figure 7

Seismicity to the east of Opak fault based on the hypocenter relocation result using jointly the data from the DOMERAPI and BMKG networks. The 2006 earthquake location was taken from global CMT catalogue (Ekström et al. 2016)

Full size image

Figure 8 shows a vertical cross section B-B′ depicting hypocenters before and after hypocenter relocation. It is shown that the dip of the slab is quite steep, similar with the one, e.g., beneath Andaman (Pesicek et al. 2008) associated with the subduction of an old oceanic lithosphere. A striking feature observed after hypocenter relocation is the nearly vertical of events down to a depth of ~30 km. We interpret that this cluster of events represents a backthrust existing to the south of central Java. A backthrust, in which its displacement is generally in an opposite direction to that of the main thrust direction, is assumed to form as a result of layer-parallel shortening in a late stage of thrust sequences (Xu et al. 2015). In the case of Java, it is similar to that existing along Sumatra, which is more obvious (Nugraha et al. 2015). This observation agrees well with the result of teletomoDD hypocenter relocation using the BMKG catalogue from 2009 to 2014 (Nugraha et al. 2015). Our interpretation is also supported by the focal mechanism of some selected events (Ekström et al. 2016) shown in Fig. 9. In general, they have a mechanism of a thrust fault that is related to the main thrust and a backthrust.

Fig. 8
figure 8

The same as Fig. 6, but for vertical cross section B-B′

Full size image
Fig. 9
figure 9

Focal mechanism solutions obtained by events selected from the global CMT catalogue (Ekström et al. 2016) showing that most events have thrust mechanisms

Full size image

4 Conclusions and further work

The hypocenter relocation result reveals some important tectonic features. One of them is the observed nearly vertical cluster of events to the south of central Java, which is related to a backthrust fault. In addition, the relocated events also depict a cluster of earthquakes parallel to the Opak fault characterized by a contrast in wavespeed (Zulfakriza et al. 2014), and a possible double seismic zone in the study region.

For further work, we will use our hypocenter relocation result to perform seismic travel time tomography in order to construct 3-D v P, v S and v P/v S models. These models will be used to understand the deep structure of the Merapi volcano, which has not been detected clearly by previous works due to lack of data coverage. As previous volcanic studies (e.g., Lei et al. 2013), a multiscale seismic tomographic imaging will be carried out utilizing the DOMERAPI and BMKG data and by incorporating the data from MERAMEX as used in previous works (Wagner et al. 2007; Koulakov et al. 2007; Rohadi et al. 2013; Haberland et al. 2014) to reconstruct magma bodies of Merapi and ascending low velocity features from the subducting slab to the volcanic arc and some other tectonic features. We will also combine data from the Center for Investigation and Technology Development of Geological Disaster, Indonesia (BPPTKG), to obtain images, especially below the Merapi summit with unprecedented detailed resolution.