ENHANCING HYPOCENTER RELOCATION FOR SEISMIC HAZARD EVALUATION IN CENTRAL JAVA AND YOGYAKARTA INDONESIA: FINDINGS FROM THE BMKG EARTHQUAKE CATALOG OF AUGUST-OCTOBER 2024

1. Introduction

Java Island is one of the regions in Indonesia with a high seismic hazard level due to its position within an active seismotectonic arc system extending from Sumatra to East Java (Aribowo et al., 2022; Muttaqy et al., 2023; Raharja et al., 2024; Soehaimi, 2008). This seismic zone is influenced by the complex interaction between the Indo Australian and Eurasian tectonic plates, with the active subduction zone south of Java serving as the primary principal of seismic activity (Gunawan and Widiyantoro, 2019; Koulali et al., 2017). Java's catastrophe index and seismic risk level are influenced by various factors, including earthquake magnitude and proximity to the seismic source, population density, local geological characteristics, and the resilience of infrastructure (Liu et al., 2024; Soehaimi, 2008). Central Java and Yogyakarta are particularly known for their relatively high levels of seismic activity, both offshore and onshore (Muttaqy et al., 2023; Ramdhan et al., 2020). Most earthquakes in this region occur in the offshore area south of Java, near the active subduction zone, indicating that tectonic plate subduction is the dominant mechanism. Conversely, shallow inland earthquakes are generally associated with the activity of local active faults. This situation highlights the importance of seismicity analysis to better understand subsurface structures, tectonic configurations, and earthquake hazard assessments.

The coexistence of subduction driven offshore earthquakes, active fault systems, and complex backarc structures creates a multihazard environment where the potential impacts extend beyond strong ground shaking to include secondary hazards such as tsunamis, liquefaction, and landslides (Aribowo et al., 2022). Thick volcanic deposits in the area can amplify seismic waves, further exacerbating structural damage in densely populated zones. These factors highlight the importance of integrating seismotectonic analysis into disaster risk reduction frameworks, ensuring that urban planning, infrastructure design, and community preparedness measures are aligned with the region’s multisource and multimechanism earthquake characteristics (Soehaimi, 2008). A comprehensive understanding of these hazards is essential for reducing potential casualties, safeguarding critical infrastructure, and enhancing the resilience of communities in Central Java and Yogyakarta. Due to the complicated interactions between tectonic processes in Yogyakarta and Central Java, thorough seismic characterization is an essential part of hazard assessment. High precision earthquake source parameters and a comprehensive seismotectonic framework are both necessary to comprehend the temporal and spatial patterns of earthquake activity. Refinement of earthquake locations by hypocenter relocation is an initial stage in this process, allowing the identification of active fault segments, elucidating subsurface structural configurations, and ultimately improving the accuracy of seismic hazard models for the area (Hidayat et al., 2024).

Accurate hypocenter location data is essential for such analyses (Hidayat et al., 2024; Muttaqy et al., 2023; Syahbana et al., 2021; Zulfakriza et al., 2024). One technique to improve the precision of earthquake localization is hypocenter relocation, which involves recalculating earthquake positions based on seismic wave travel time data, thereby reducing errors related to the velocity model used (Garini et al., 2014). Accurate hypocenter relocation enables better visualization of earthquake distributions, which can reflect underlying geological structures and provide deeper insight into regional seismotectonic environments (Waldhauser and Ellsworth, 2000). Hypocenter relocation is a technique for recalculating or correcting the position of hypocenter using seismic data to achieve a more precise location (Garini et al., 2014). Since it enables the viewing of hypocenter continuity that reflects subsurface structures, precise hypocenter relocation is essential for increasing earthquake localization precision (Waldhauser, 2001). Additionally, it makes it possible to identify seismotectonic settings by looking at the spatial distribution of earthquakes. The Double-Difference (DD) method is one of the most used approaches for hypocenter displacement. The relative locations of hypocenters are ascertained using the DD approach (Hidayat et al., 2022, 2024; Waldhauser, 2001). It is based on the variations in travel times between earthquake pairs that were recorded at the same seismic station. The principle states that if the distance between two paired earthquakes is relatively modest in comparison to their distances from the seismic station, their ray trajectories and waveforms can be deemed substantially equal. This assumption minimizes mistakes caused by the velocity model by allowing the travel time difference between the two earthquakes recorded at the same station to be interpreted as a function of the distance between their hypocenters (Waldhauser dan Ellsworth, 2000).

Central Java and Yogyakarta have been identified as the most densely inhabited regions in Indonesia, with over 73 million inhabitants residing in this highly seismic zone (Muttaqy et al., 2023). Given the significant seismic risk, the examination of earthquake clusters in this area is crucial to enhance and highlight the Indonesian seismic hazard map (Muttaqy et al., 2023). We performed hypocenter relocation using the 1-D AK135 velocity model, following the methodology previously applied in studies focusing on the Central Java and Yogyakarta regions (Ramdhan et al., 2020; Saputra et al., 2021; Wagner et al., 2007). The primary dataset for this relocation process was the earthquake event catalog data from August to October 2024, which we obtained from the BMKG Banjarnegara Station and processed using SeisComp4. This research aims to analyze the seismicity of Central Java and Yogyakarta in both spatial and temporal dimensions, with particular emphasis on the subduction zone and its associated structural features. In this research, we considered the existence of active earthquake sources in the Yogyakarta and Central Java regions that have not been fully identified by the National Center for Earthquake Studies (PUSGEN). By considering this, it is expected that this research can produce new findings related to the distribution of local earthquake activity and indications of active fault structures that have not been documented before, thereby contributing to updating disaster information in the region. The results are expected to provide a foundational reference for future, more comprehensive studies and to serve as an initial step toward effective earthquake disaster mitigation in the region.

Refining earthquake hypocenter locations is not just a methodological improvement; it is a basic need to make seismic source characterization less uncertain. Accurate relocation helps to map out active fault segments, make the geometry of subduction-related seismicity clearer, and find hidden seismotectonic structures that normal catalog data can't show. These improvements directly help make seismic hazard models more accurate in places like Central Java and Yogyakarta, which are very densely populated and culturally important. This makes mitigation strategies stronger, helps with urban planning, and makes communities more resilient to future earthquakes.

2. Geodinamics of Central Java and Yogyakarta

The geodynamics of Central Java and Yogyakarta are influenced by the interplay of subduction processes to the south of Java, active volcanism to the north, and fault systems that extend E – W and SW – NE (see Figure 1). Because of how the tectonic plates are arranged, the area is more likely to have earthquakes, both shallow and intermediate, that can do a lot of damage in locations with a lot of people. The presence of unconsolidated volcanic material, basin sedimentary deposits, and steep terrain further exacerbates the seismic hazard, raising ground shaking amplification and increasing the possibility for liquefaction in specific zones. These conditions illustrate that regional geodynamics not only modify the terrain and geological history, but they also make earthquakes more likely, which can have a huge influence on infrastructure and society in this area.

The Indo-Australian oceanic plate is moving south of Java and sliding beneath the Eurasian continental plate. This is causing the area around Central Java and Yogyakarta to change over time. This ongoing subduction has created the Central Java Quaternary Volcanic Arc, which is located to several active volcanoes, such as Merapi, Merbabu, Sindoro, Sumbing, and Slamet. These volcanoes are still showing signs of magmatic activity. The Southern Mountains, made up of mostly Tertiary carbonate rocks, are the most prominent feature in the southern part of the region. They make up the famous Gunung Sewu karst system, which has many caves and rivers that run underground. The Kendeng Zone to the north is made up of sedimentary rocks from the Miocene to the Plio-Pleistocene that have been bent into complicated folds and faults. The Opak Fault and other active fault systems show that tectonic processes are just as important as volcanism in shaping the landscape and hazards in the area. All of these things together make Central Java and Yogyakarta two of the most seismically active places in Indonesia. The area has tectonic and volcanic earthquakes that happen over and over again. One of the most damaging was the 2006 Yogyakarta earthquake. Unconsolidated volcanic deposits and sedimentary basins make the area even more dangerous because they tend to make the ground shake more and cause secondary hazards like liquefaction and slope failures.

Figure 1. Geological map of the Central Java sheet at a scale of 1:500.000 (modified from Amin et al., (1999))

2. Methods

We used the Double Difference (DD) method for this study. The DD approach is highly effective in areas with numerous closely situated local earthquake events clusters, resulting in a more concentrated distribution that aligns with the geological structure (Muttaqy et al., 2023; Zhang, 2003; Zhang and Thurber, 2006). The DD method relies on the assumption that when the distance between two earthquake sources is considerably smaller than their distances to a recording station, and the velocity heterogeneity along their respective ray paths is roughly equivalent, the wave rays from both sources can be regarded as nearly indistinguishable (Waldhauser and Ellsworth, 2000). The precision of hypocenter determination is affected by various parameters, such as the configuration of the seismic station network, the quality and quantity of seismic phases, and the accuracy of the recorded arrival timings. The DD technique posits that when the interhypocentral distance of two earthquakes is minimal compared to their distance from the station, their wave raypaths reaching the station are substantially identical. Consequently, the disparity in transit time between the two events, measured at a shared station, can be understood as a function of the relative distance between their hypocenters (Hidayat et al., 2022). This method enables the minimization of velocity model error without requiring station corrections (Waldhauser and Ellsworth, 2000). The DD method has certain limitations, particularly when earthquake distributions are sparse or do not form clusters, making the principle of comparing nearby earthquake pairs less effective. Moreover, DD relocation requires handling a large number of earthquake pairs, which results in higher computational demands, especially for large datasets (Hidayat et al., 2024).

The travel time of seismic waves from an earthquake source iii to a seismic station kkk can be expressed using ray theory (Um & Thurber, 1987):

[TeX:] T_{k}^{i} = \tau^{i} + \int_{i}^{k}{u\ ds}$T_k^i={\tau }^i+{\int }_i^{k}uds$ (1)

where [TeX:] \tau^{i}${\tau }^i$ represents the origin time of the earthquake, [TeX:] i,\ u$i,u$ is the slowness vector, and [TeX:] ds$ds$ denotes the partition of the ray path. The relationship between arrival time and source location is highly non-linear; therefore, linearization is applied by calculating the misfit between observed and predicted travel times, given by:

[TeX:] r_{k}^{i} = \sum_{l = 1}^{3}{\frac{\partial T_{k}^{i}}{\partial x_{l}^{i}}\mathrm{\Delta}x_{l}^{i}} + \mathrm{\Delta}\tau^{i} + \int_{i}^{k}{\delta u\ ds}$r_k^i={\sum }_{l=1}^3\frac{\partial T_k^i}{\partial x_l^i}\mathrm{\Delta }{x}_l^i+\mathrm{\Delta }{\tau }^i+{\int }_i^k\delta uds$ (2)

By substituting the arrival times of another earthquake source [TeX:] j$j$ recorded at the same station [TeX:] k$k$, the expression becomes:

[TeX:] r_{k}^{i} - r_{k}^{j} = \sum_{l = 1}^{3}{\frac{\partial T_{k}^{i}}{\partial x_{l}^{i}}\mathrm{\Delta}x_{l}^{i}} + \mathrm{\Delta}\tau^{i} + \int_{i}^{k}{\delta u\ ds} - \sum_{l = 1}^{3}{\frac{\partial T_{k}^{j}}{\partial x_{l}^{j}}\mathrm{\Delta}x_{l}^{j}} - \mathrm{\Delta}\tau^{j}$r_k^i-r_k^j={\sum }_{l=1}^3\frac{\partial T_k^i}{\partial x_l^i}\mathrm{\Delta }{x}_l^i+\mathrm{\Delta }{\tau }^i+{\int }_i^k\delta uds-{\sum }_{l=1}^3\frac{\partial T_k^j}{\partial x_l^j}\mathrm{\Delta }{x}_l^j-\mathrm{\Delta }{\tau }^j$

[TeX:] - \int_{i}^{j}{\delta u\ ds}$-{\int }_i^j\delta uds$. (3)

[TeX:] r_{k}^{i} - r_{k}^{j} = \sum_{l = 1}^{3}{\frac{\partial T_{k}^{i}}{\partial x_{l}^{i}}\mathrm{\Delta}x_{l}^{i}} + \mathrm{\Delta}\tau^{i} - \sum_{l = 1}^{3}{\frac{\partial T_{k}^{j}}{\partial x_{l}^{j}}\mathrm{\Delta}x_{l}^{j}} - \mathrm{\Delta}\tau^{j}$r_k^i-r_k^j={\sum }_{l=1}^3\frac{\partial T_k^i}{\partial x_l^i}\mathrm{\Delta }{x}_l^i+\mathrm{\Delta }{\tau }^i-{\sum }_{l=1}^3\frac{\partial T_k^j}{\partial x_l^j}\mathrm{\Delta }{x}_l^j-\mathrm{\Delta }{\tau }^j$ (4)

By assuming that the two sources are located very close to each other, the ray paths become nearly identical and cancel out, reducing Equation 3 to Equation 4, where [TeX:] r_{k}^{i} - r_{k}^{j}$r_k^i-r_k^j$ represents the Double-Difference (Waldhauser and Ellsworth, 2000). The Double-Difference method refers to the difference in travel times between observed data and calculated data for two nearby earthquake sources, which can be expressed as:

[TeX:] r_{k}^{i} - r_{k}^{j} = {{(T}_{k}^{i} - T_{k}^{j})}^{obs} - {{(T}_{k}^{i} - T_{k}^{j})}^{cal}$r_k^i-r_k^j={{(T}_k^i-T_k^j)}^{obs}-{{(T}_k^i-T_k^j)}^{cal}$ (5)

To further investigate the spatial orientation of these relocated hypocenters, it is necessary to compute the azimuth and relative distance between earthquake pairs. This is achieved using the Haversine formula, which provides an accurate estimation of great-circle distances on a spherical Earth:

[TeX:] d = 2R\ arcsin\left( \sqrt{\sin^{2}\left( \frac{\mathrm{\Delta}\varphi}{2} \right) + cos\left( \varphi_{1} \right)\ cos\left( \varphi_{2} \right){\ sin}^{2}\left( \frac{\mathrm{\Delta}\lambda}{2} \right)} \right)$d=2Rarcsin\left(\sqrt{{sin}^2\left(\frac{\mathrm{\Delta }\phi }{2}\right)+cos\left({\phi }_1\right)cos\left({\phi }_2\right){sin}^2\left(\frac{\mathrm{\Delta }\lambda }{2}\right)}\right)$ (6)

where [TeX:] R$R$ is the Earth’s radius, [TeX:] \Delta\varphi = \varphi_{2} - \ \varphi_{1}$\Delta \phi ={\phi }_2-{\phi }_1$ is the difference in latitude, and [TeX:] \Delta\varphi = \lambda_{2} - \ \lambda_{1}$\Delta \phi ={\lambda }_2-{\lambda }_1$ is the difference in longitude between two hypocenters. The calculated azimuthal orientations are then represented using a rose diagram, which provides a statistical visualization of dominant structural trends, seismic alignments, and fault orientations. In this study, the integration of double-difference relocation with Haversine based azimuth analysis allows for a more robust characterization of seismicity patterns in Central Java and Yogyakarta.

The data used for this research came from a catalog of earthquakes in Central Java and Yogyakarta from August to October 2024 and from stations in Banjarnegara using Seiscomp-4 with a geographic range of 6.0˚-11.0˚ and 108.5˚-111.6˚ (see Figure 2). The data catalog comprises information about the earthquake occurrence, date, time, location, depth, agency, mode, status, RMS residual, azimuthal gap, magnitude, phase, and the station that recorded the earthquake.

Figure 2. The distribution of the seismic network in Central Java and Yogyakarta is represented by blue triangle symbols, while the red triangles indicate volcanoes located within the study area. All elements are overlaid with the active fault map in the region

The hypocenter relocated in this research through two programs, ph2dt and hypoDD. As required by the relocation software, the first step was to convert the original BMKG earthquake catalog data from.txt format to.pha format. The (*.pha) catalog file is the primary input for the ph2dt program. Important details about the event, including the earthquake's date and time, hypocenter coordinates (latitude and longitude), depth, magnitude, RMS error, seismic station, seismic wave arrival times, quality indicators, and wave phases, are provided by the (*.pha) catalog. We calculate the difference in travel times between earthquake pairs using this precise event information. The hypocenter is then moved using the hypoDD algorithm. The hypoDD program's Double-Difference (DD) method for moving things around entails putting similar earthquake occurrences into clusters based on a set distance threshold. The MINLNK parameter can be used to find the best arrangement of these clusters by setting the least number of observations needed to form a cluster equal to the minimum number of observations needed for events that are nearest through. The initial development of clusters is additionally affected by the maximum inter-event distance within a cluster, denoted by the WDCT parameter. To ensure dependable outcomes, the WDCT value must be equal to or less than the MAXSEP parameter, which delineates the maximum permissible spacing between event pairings. A further crucial factor is the implementation of a suitable one-dimensional (1D) velocity model. This investigation utilized the AK135 velocity model, characterized by a Vp/Vs ratio of 1.74, as illustrated in Figure 2. This model establishes the foundation for computing travel durations and is essential for the precision of relocations. Furthermore, an essential parameter that requires meticulous assessment is the Condition Number (CND), which indicates the quality of the inversion process. For dependable relocation outcomes, CND values must reside between 40 and 80. A low CND value may signify overdamping, wherein excessive regularization causes few or no alterations in event locations throughout the inversion process, leading to trivial decreases in travel-time residuals. A high CND number may indicate underdamping, where insufficient regularization leads to unstable and excessive shifts in event positions, thereby jeopardizing the accuracy of the relocation results.

Figure 3. Velocity profiles of P-waves (Vp) and S-waves (Vs) based on the AK135 velocity model down to a depth of 700 km. The blue line represents S-wave velocities, whereas the red line represents P-wave velocities

The dataset used in this study consists of earthquake catalog data from BMKG Banjarnegara, covering the Central Java and Yogyakarta regions for the period of August to October 2024. A total of 393 earthquake events were recorded in August 2024, 424 events in September 2024, and 276 events in October 2024, resulting in a total of 1,093 events recorded by 29 seismic stations. The recorded earthquake magnitudes ranged from 1.0 to 5.8 Mw, with focal depths ranging from 5 to 256 km. The earthquake data were initially processed using SeisComP4 and then relocated using the Double-Difference (DD) method. The relocation process involved two main steps: using the ph2dt program followed by the hypoDD program.

Relocation process involved two main steps: first, applying the ph2dt program to calculate differential travel times, and then using the hypoDD program to relocate the earthquake hypocenters. The number of earthquake events before and after relocation differed because not all events satisfied the selection criteria required for processing in ph2dt. Specifically, events were excluded if they had an insufficient number of associated phase readings, poor azimuthal station coverage, or large residuals in the travel-time data. These conditions are necessary to ensure that the double-difference algorithm can accurately minimize relative location errors between nearby event pairs. In addition, some relocated events exhibited so-called “airquake” effects, in which the hypocenter positions were incorrectly shifted to locations above the Earth’s surface. Such artifacts typically occur when the available phase data are too sparse or unevenly distributed, leading to unstable inversion results in hypoDD. Therefore, only events that fulfilled the minimum criteria of station coverage, phase quality, and residual thresholds were retained for the final relocated catalog.

Hypocenter relocation in this study was carried out using HypoDD. The parameters applied in the hypocenter relocation process for August, September, and October 2024 included MAXNGH, MINWGHT, MINLNK, MINOBS, and MAXOBS. The MAXNGH parameter was consistently set to 16 for each month. The maximum distance between events and stations (MAXDIST) was set at 600 km, while MAXSEP, which represents the maximum separation between event pairs, was fixed at 100 km throughout the three-month observation period. The selection of these parameters was meant to make sure that the inversion process was stable and that the relocated hypocenters were reliable. We chose a MAXNGH value of 16 to keep enough neighbouring events in each cluster without putting too much strain on the computer. The relatively large MAXDIST (600 km) made it possible to include stations that were far away and gave more phase data, which improved azimuthal coverage and made it less likely that the location was wrong. The MAXSEP of 100 km made sure that event pairs were still related in space, but it also let bigger clusters form in areas with few earthquakes. Similarly, thresholds for MINWGHT, MINLNK, MINOBS, and MAXOBS were applied to guarantee that only events with adequate observational constraints and reasonable travel-time residuals were included in the relocation process. These settings reflect a balance between data quality control and the need to retain a sufficient number of events for meaningful seismicity analysis. The ph2dt program provides output in the form of dt.ct data containing the travel time differences of the selected event pairs, event.sel containing information on the events selected as pairs, and ph2dt.log containing information on the results of the ph2dt process. The output of the ph2dt program is used for the relocation process in the hypoDD program. There are important things in the ph2dt program, namely outliers and weakly linked events. The program will reduce events that are considered outliers, the smaller the percentage of outliers, the better the data. Furthermore, only events that have a strong link will be used in the relocation process, so the smaller the percentage of weakly linked events, the fewer weakly linked event pairs there will be.

The proportion of data classified as outliers was relatively low, accounting for 1% in August (756 pairs) and September (843 pairs), and 0% in October (266 pairs). After filtering, the number of P-wave phase pairs used for relocation was 13,801 (54%) in August, 15,100 (54%) in September, and 9,968 (55%) in October. Likewise, the number of S-wave phase pairs used was 12,301 (56%), 11,230 (54%), and 8,835 (57%), respectively. Regarding the connectivity between events, a number of earthquakes were found to be weakly connected, with 33 events (8%) in August, 21 events (4%) in September, and 17 events (6%) in October. The total number of earthquake pairs formed was 4,114 (August), 4,497 (September), and 2,938 (October), with an average number of phase records per pair being 6, 5, and 6, respectively.

Figure 4 depicts the distribution of earthquake hypocenter positions in August 2024, September 2024, and October 2024. The green circles on the maps reflect the hypocenter positions before and after relocation, whereas the yellow circles represent the hypocenter positions following relocation. Each circle's size corresponds to the magnitude of the earthquake. Magnitudes range from 1 Mw to ≥ 4 Mw, as represented by varying circle sizes. Volcanic locations are indicated by red triangles, whereas seismic stations are represented by blue triangles. Figure 4a illustrates the distribution of seismic events in August 2024, Figure 4b depicts the distribution for September 2024, and Figure 4c represents the distribution for October 2024. In August 2024, 393 earthquake events were initially recorded. After processing with ph2dt, the dataset was reduced to 392 events, and after relocation with hypoDD, 354 events remained. In September 2024, 424 events were initially recorded, which became 421 after ph2dt processing and 391 after hypoDD relocation. In October 2024, the original 276 events were reduced to 275 after ph2dt and finally to 257 events after relocation using hypoDD.

Figure 4. Spatial distribution of earthquake events during (a) August 2024, (b) September 2024, and (c) October 2024. Earthquake depth is colour-coded: yellow for shallow, green for moderate, and orange for deep events. The blue triangles indicate the locations of seismic stations, and the red triangles represent volcanoes.

3. Results and Discussion

3.1. Distribution of Hypocenter Locations Before and After Relocation

The distribution of hypocenters during the observation period is shown in the Comparison Maps of Hypocenter Positions Before and After Relocation for August, September, and October 2024 (see Figure 5). In these maps, the green circles illustrate the hypocenter locations before relocation, whereas the yellow circles denote the positions following the relocation process. The diameter of each circle reflects the magnitude of the earthquake, with larger circles indicating magnitudes of 1 Mw, 2 Mw, 3 Mw, and ≥ 4 Mw, respectively. In the study area, the red triangular symbols mark the active volcanoes, whereas the blue triangles represent the seismic stations utilized for data collection. This visualization effectively showcases the spatial distribution and magnitude of seismic events, while also delivering crucial insight for analyzing earthquake occurrence patterns and their connection to geological structures. Additionally, it highlights the distribution of seismic stations, which plays a vital role in enhancing the accuracy of hypocenter determination. The map illustrates the hypocenter's position before and after the relocation, suggesting that the hypocenter seems more concentrated and centralized. A rose diagram, as seen in Figure 6, is constructed to determine the direction of the shift in hypocenter relocation.

Figure 5. Spatial distribution map of hypocenters before (green circles) and after (yellow circles) relocation. The size of each circle is proportional to the moment magnitude of the event.

The depiction of the shift pattern of the earthquake hypocenter position after being relocated by the DD method can be seen in Figure 6. The data used in making this rose diagram are latitude and longitude data before relocation with latitude and longitude data after relocation. The distance and azimuth shift values can be calculated using the Haversine formula. Figure 6 shows the shift in hypocenter position after relocation in August 2024, September 2024, and October 2024. Green to yellow colours indicate the distance of the hypocenter position shift after relocation, from 0 km to > 60 km. The diameter of the circle indicates the frequency of the data in its azimuthal direction (N, N-E, E, S-E, S, S-W, W, and N-W). Figure 6 shows the rose diagram of the shift of the hypocenter position after relocation. The rose diagram shows that the position of the hypocenter has shifted in all directions. In August 2024, the hypocenter position shift tends to azimuth between 85˚ - 95˚ to the east and azimuth between 265˚ - 275˚ to the west with a shift distance of 0 to > 41 km. In September 2024, the hypocenter position shift tends to azimuth between 50˚ - 104˚ to the east and azimuth between 240˚ - 280˚ to the west with a shift distance of 0 to > 49 km. In October 2024, the hypocenter position shift tends to azimuth between 50˚ - 110˚ to the east and azimuth between 265˚ - 275˚ to the west with a shift distance of 0 to > 60 km. So that the tendency of shifting the position of the hypocenter for the period August-October 2024 tends to shift to the E- W direction

Figure 6. Rose diagrams of hypocenter location shifting after relocation, with colour indicating the displacement distance: (a) August 2024, (b) September 2024, and (c) October 2024.

Figure 6 illustrates a comparative analysis of travel time residuals before and after the relocation process. The construction of the residual histogram aims to evaluate the accuracy of the relocation results. A residual value that approaches zero indicates that the calculated travel time is close to the actual travel time, which in turn suggests that the relocated hypocenter is closer to its true position. Greyscale bars represent the residuals before relocation, while black bars indicate the residuals after relocation. Each month exhibits different residual frequency distributions, reflecting variations in data quality and event characteristics over time. The data used in generating these histograms consist of travel time residuals both before and after the relocation process. Figure 7a presents the travel-time residuals of the hypocenter before and after relocation for August 2024. The histogram shows that the residual values prior to relocation ranged from -3 s to 3 s. The frequency of residuals between -1 s and 1 s was approximately 18,000, while the frequency of residuals less than -1 s and greater than 1 s was around 3,800. After relocation, the residual values ranged from -1.7 s to 1.7 s. The frequency of residuals between -1s and 1s was approximately 16,000, whereas those beyond this range totaled around 1,900. Figure 7b shows the travel-time residuals of earthquakes before and after relocation for September 2024. The histogram indicates that the residuals before relocation ranged from -4 s to 4 s. The frequency of residuals within -1 s to 1 s was about 18,000, while those outside this range were approximately 3,800. After relocation, the residual values were within -1.9 s to 1.9 s, with around 18,000 residuals falling within -1 s to 1 s and only about 900 beyond that range. Figure 7c illustrates the travel-time residuals before and after relocation for October 2024. Before relocation, the residuals ranged from -4 s to 4 s. The number of residuals between -1 s and 1 s was approximately 13,000, while those outside this range numbered around 3,800. After relocation, the residual range narrowed to -2 s to 2 s, with around 16,000 residuals between -1 s and 1 s and about 2,500 outside that range. The results after relocation show better accuracy for the earthquake hypocenter, as the data shows that the residuals are more closely grouped around zero in every month. This implies that the positions of the earthquake hypocenters are getting closer to their actual locations (Hidayat et al., 2021).

Figure 7. Histogram Comparison of Travel Time Residuals Before and After Relocation: (a) RMS August 2024, (b) RMS September 2024, and (C) RMS October 2024. Greyscale bars represent the travel time residuals before relocation, while black bars indicate the travel time residuals after relocation.

3.3. Seismicity Characteristics of the Central Java and Yogyakarta Region

Figure 8 displays the relocated hypocenter locations. The size of the circular symbols represents the earthquake magnitude, while the colour indicates the depth. From the smallest to the largest, the circles correspond to magnitudes of 1 Mw, 2 Mw, 3 Mw, and ≥ 4 Mw, respectively. Yellow indicates shallow earthquakes (0-20 km), green represents moderate depth earthquakes (20-60 km), and orange indicates deep earthquakes (> 60 km). The red triangles mark the locations of volcanoes, while the blue triangles represent the positions of seismic stations. We observed fluctuations in the number of earthquake events over the course of three months in the Central Java and Yogyakarta regions. In August, a total of 354 earthquake events were recorded, consisting of 121 shallow earthquakes, 189 moderate depth earthquakes, and 44 deep earthquakes. In September, there was an increase in the number of events, particularly in shallow earthquakes, which rose from 121 in August to 161 events. In contrast, the number of moderate depth earthquakes decreased from 189 to 175 events, while deep earthquakes increased from 44 to 55 events. The total number of earthquake events in September reached 391. In October, a significant decrease in the number of events was observed, especially in shallow and moderate depth events. A total of 257 earthquakes were recorded during this month, comprising 101 shallow events, 103 moderate depth events, and 53 deep events. The spatial distribution of earthquakes is primarily concentrated in the Central Java and Yogyakarta regions, especially along the Baribis - Kendeng Fault (yellow box) and the Opak Fault Zone (blue box), according to active fault data from the National Center for Earthquake Studies (Damanik et al., 2021; Gunawan and Widiyantoro, 2019; Novianto et al., 2021; Pratama et al., 2021; Saputra et al., 2021; Tsuji et al., 2009). Furthermore, there are clusters of seismic occurrences that have not been correlated with any recognized active fault, indicated by the red boxes (see Figure 8).

Figure 8. Seismicity maps based on relocated hypocenters in the Central Java and Yogyakarta regions are presented as follows: a) seismicity in August, b) seismicity in September, and c) seismicity in October 2024. The seismicity maps are overlaid with a cross-section A-A’. Yellow dots represent shallow earthquakes, green indicates moderate depth events, and orange corresponds to deep earthquakes. Blue triangles denote the distribution of seismic stations, while red triangles indicate the locations of volcanoes.

Seismic activity is prominently detected in the Baribis - Kendeng Fault zone. The existence of an active Baribis Fault extending eastward from the Sunda Strait together with a Kendeng Fault traversing from West Java to the east (Damanik et al., 2021; Gunawan and Widiyantoro, 2019). This fault extends from East Java to West Java, the Kendeng Fault segment dominates in the eastern and central parts, while the Baribis Fault is in the western part. The Kendeng-Baribis Fault is part of a back-arc thrust system in northern Java formed by the collision of the Indo-Australian Plate with the Eurasian Plate (Koulali et al., 2017; Raharja et al., 2024). The results of a recent study conducted by Damanik et al., (2021), mainly using a network of borehole seismometers deployed around the Baribis segment from 2019 to 2021, have identified a number of small earthquake events with magnitudes of 2 - 4 Mw that are closely located, showing an earthquake event on December 10, 2019, with a magnitude moment of 3.0 Mw and on March 29, 2020, of 3.2 Mw. An earthquake occurred on August 3, 2024, with a moment magnitude of 2.3 Mw, situated north of the Baribis - Kendeng Fault (see Figure 8a and Figure 8c). Subsequently, a comparable event transpired south of the same fault zone on October 20, 2024, registering a moment magnitude of 2.6 Mw. These results indicate persistent tectonic activity along the Baribis - Kendeng Fault. Despite the relatively small magnitude, the event's location and mechanism furnish direct evidence of the seismotectonic activity of the Baribis - Kendeng Fault segment.

The upcoming seismic cluster is located in the Central Java and Yogyakarta regions, with a SW-NE orientation. Administratively, it covers Kebumen Central Java and Kulon Progo Yogyakarta. Based on the direction and distribution pattern of earthquake occurrences, we assume that the cluster in the Kebumen area is associated with the Meratus Pattern, also known as the Kebumen - Meratus Trend (Satyana, 2008), which extends as far as Kalimantan. Other researchers have referred to this fault as the Central Java Fault (CJF) (Hall et al., 2007; Hoffmann-Rothe et al., 2001). Seismic activity in this area is relatively intense, with earthquakes occurring at shallow to moderate depths. Satyana (2008) identified the Muria-Kebumen Fault in Central Java as one of the main Paleogene structures that has a SW-NE orientation and left-lateral strike-slip characteristics. This structure is also referred to as the Meratus Trend due to its alignment and possible tectonic linkage with the structural line of the Meratus Mountains in Kalimantan. From a seismological perspective, this SW-NE orientation has the potential to play a role in controlling the direction of propagation and segmentation of earthquake sources in the southern region of Central Java. The tectonic linkage between the Meratus Trend indicates a structural continuity formed since the Paleogene, which then experienced reactivation due to Quaternary tectonic forces. Seismic activity was recorded in August 2024 and continued to decline until October 2024, with a moment magnitude > 3 Mw with shallow to moderate depth. The results of research from Novianto et al. (2025) showed a low Bouguer anomaly value in the direction of the Kebumen - Meratus Trend. The evidence of this earthquake cluster reinforces that this zone is likely to contain active faults that have not been mapped by the National Earthquake Study Center.

Moving eastward, specifically in the Kulon Progo region of Yogyakarta, there is a noticeable concentration of earthquake events with a dominant SW - NE orientation. These events were particularly active during August and September 2025. Unlike the cluster located further west, this cluster consists predominantly of deep focus earthquakes. Previous studies have identified this cluster as part of the Progo - Muria Lineament adapted from (Hall et al., 2007; Novianto et al., 2025), a tectonic lineament that extends from south to north and continues toward the Muria Mountains. This lineament is referred to as the inferred Progo - Muria Lineament or the strikeslip Central Java Fault. The frequency of seismic occurrences in the Kulon Progo area of Yogyakarta signifies a busy tectonic zone characterized by intricate geodynamics. Deep earthquakes dominate this cluster, contrasting with shallow seismicity in the surrounding area, and indicate the involvement of subcrustal processes that may relate to the geometry of the subduction slab beneath southern Central Java. The temporal clustering of seismic events in August and September 2025 confirms the potential for episodic stress release patterns, which may be influenced by regional plate convergence and local fault kinematics.

From a hazard perspective, the Progo - Muria Linement has the potential to generate medium to large magnitude earthquakes with a significant horizontal movement component, which can create building damage. Although deep-focused earthquakes generally produce lower shaking intensity at the surface than shallow earthquakes, their proximity to densely populated areas increases the potential for cumulative hazards, especially if accompanied by possible triggering of shallow seismicity. The SW-NE orientation of this seismic cluster corresponds with the structural directions of several active faults in Central Java, enhancing concerns regarding the potential for fault to fault rupture propagation during significant events. Due to its tectonic significance and recurrent seismic activity, the Progo - Muria Linement line should be considered one of the main sources of seismic hazards in regional disaster risk assessments. This requires continuous seismic monitoring to identify changes in activity rate or hypocenter depth distribution, detailed fault mapping and paleoseismological studies to estimate its return interval and shear rate, and integration of this fault line into a probabilistic seismic hazard model for Central Java.

The persistent seismic pattern east of the Opak Fault Zone, characterized by earthquake events distributed approximately 5-10 km from its primary trace and aligned SW-NE, demonstrates the intricate nature of fault interactions in this area. While Ramdhan et al. (2020) noted that the Opak Fault Zone inclines to the east, cross-sectional analysis indicates that the recorded seismicity does not exhibit a distinct cluster directly on the fault plane. Instead, the eastward distribution of hypocenters points to the possible involvement of secondary fault segments or other structures in the vicinity, such as the Oyo Fault or Ngalang Fault (see Figure 9). From a hazard perspective, the proximity of these seismic sources to the Opak Fault Zone indicates the possibility of an interconnected fault network with the potential to produce damaging earthquakes. Multiple fault interactions, including stress transfer from one segment to another, may increase the probability of complex rifting scenarios extending beyond the main line of the Opak Fault Zone. The circumstances could contribute to increased fault lengths, higher magnitudes, and extended areas of strong shaking, particularly if the displacement crosses the fault border.

The location of seismic events in areas with high population density has direct implications for increasing seismic risk levels. Residential areas located along the Opak Fault Zone covering the Bantul, Sleman, and Gunungkidul regencies have a significant level of vulnerability. The presence of thick volcanic deposits around Yogyakarta has the potential to cause seismic wave amplification, which in turn can magnify the intensity of shaking and exacerbate infrastructure damage. Analysis of the events following the 2006 Yogyakarta earthquake shows that fatalities were generally caused by building collapses rather than by the earthquake shaking itself. Meanwhile, most of the injuries occurred as a result of mass panic and lack of preparedness. This confirms that structural and non-structural vulnerability factors play a major role in determining the level of damage and the number of casualties. Within the framework of disaster risk management, adequate investment is needed in Disaster Risk Reduction (DRR) programs that include improving the quality of earthquake-resistant buildings, hazard-based spatial planning, and community education. In areas that have been affected by disasters, the application of the build back better and safer principle is very important in the rehabilitation and reconstruction phases to ensure that the rebuilt infrastructure is more resilient to future seismic events.

Figure 9. Vertical cross-section along profile B-B’ depicting seismicity distribution for (a) August 2024, (b) September 2024, and (c) October 2024 inside the Opak Fault Zone section and surrounding areas.

The A-A’ cross-sectional profiles for August, September, and October 2024 are presented in Figure 7. These profiles are correlated with the Slab 1.0 model (Hayes et al., 2012). Figure 8 shows a vertical cross-section along A-A’, cutting through the subduction zone from south to north and intersecting the Opak Fault Zone. The results reveal a distribution of hypocenters that follows a sharply dipping subducting slab. This steep slab geometry is likely associated with the relatively old age of the oceanic plate (Pesicek et al., 2008). In September, a concentration of seismic events is observed at a depth of approximately 50 km, indicated by red circles. This pattern corresponds to increased seismicity in the offshore region south of Java during that month, as shown in Figure 10. The concentration of seismicity at depths of roughly 50 km, which coincides with increased activity off the southern coast of Java, indicating a risk from intermediate crustal earthquakes. Earthquakes at this depths can generate significant shaking in coastal areas, including the potential for seabed deformation and local tsunamis, especially if the process involves upward fault movement in accretionary prisms (Aribowo et al., 2022; Xu et al., 2015). The identification of shallow earthquake clusters associated with backthrust structures also adds to the complexity of the hazard source. Although generally relatively small in size, backthrust earthquakes occurring near the surface can cause high shaking intensity in the coastal zone and southern mainland of Yogyakarta, given their proximity to the population. In the Java subduction zone, the main thrust dips northward, consistent with the dominant tectonic stress orientation across the island. Between depths of 50 and 250 km, a cluster of earthquakes is distributed along the slab interface. These are intraslab earthquakes, occurring within the subducting oceanic plate along the Java Benioff Zone. Their distribution forms an inclined pattern that aligns with the subduction dip. Such seismicity is generally attributed to slab-pull forces acting within the descending plate (Lallemand et al., 2005).

Figure 10. Vertical cross-sections along profile A-A’ depicting the seismicity distribution for (a) August 2024, (b) September 2024, and (c) October 2024, overlay on the Slab 1.0 model (Hayes et al., 2012).

According to the discussion, it can be concluded that the Central Java and Yogyakarta regions are located in a complex tectonic environment, influenced by the interaction of multiple earthquake sources ranging from the main subduction zone, shallow crustal faults such as Baribis - Kendeng and Progo - Muria Lineament, to local faults around Opak Fault Zone. According to our analysis of the Progo - Muria lineament, seismic activity in August and September belonged to the category of deep earthquakes, which is different from the pattern on the Kebumen-Meratus Trend that is dominated by shallow to medium earthquakes. We identified that the earthquake events in that period were more caused by subduction/deep earthquake activity, and not from the Progo - Muria lineament. Hence, a seismic data sequence with a longer recording period is needed to confirm whether the earthquakes actually originated from the Progo - Muria lineament. The presence of shallow earthquakes on accretionary prisms, intermediate crustal earthquakes at depths of around 50 km, and intraslab earthquakes at depths >50 km creates a diverse spectrum of threats, both in terms of shaking intensity, potential infrastructure damage, and tsunami risk. Seismicity patterns distributed according to regional structural orientation show the existence of linkages between earthquake sources that have the potential to trigger compound fault scenarios. This condition requires the integration of the latest seismotectonic data, modelling of multisource shaking scenarios, and detailed mapping of active fault lines in the preparation of seismic hazard maps. In addition, disaster mitigation should focus on improving building resilience, strengthening hazard-based spatial planning, and educating the public to minimize the risk of casualties and material losses when large earthquakes occur in the future.

Subduction zone earthquakes have a significant potential to trigger liquefaction, particularly in coastal and lowland areas characterized by loose sediments, water saturated sands, and alluvial deposits. Liquefaction occurs when pore water pressure increases due to dynamic earthquake loading, causing the soil to lose its shear strength and behave like a fluid The 2018 Palu earthquake serves as a notable example in Indonesia, where widespread liquefaction induced massive lateral ground displacements and extensive infrastructure damage (Cilia et al., 2021; Susanta et al., 2019; Upomo et al., 2023). Similarly, research in Bengkulu has shown that the 2007 subduction earthquake had the potential to trigger shallow liquefaction in coastal zones (Mase et al., 2024; Mase et al., 2025). As shown by research in Mae Lao and Bengkulu, soil conditions with low shear wave velocities (Vs) have been shown to increase seismic wave amplification and liquefaction potential in shallow water-saturated sand layers (Mase et al., 2020; Mase et al., 2024). Similar phenomena are particularly relevant in Yogyakarta, which is composed of young volcanic deposits and unconsolidated alluvium, making it vulnerable to ground deformation during large earthquakes. Therefore, seismic hazard mitigation in this region needs to be conducted comprehensively, including microtremor measurements, horizontal to vertical spectral ratio (HVSR) analysis, shear wave velocity inversion, as well as microzonation mapping to determine liquefaction-prone and amplification zones. From a hydrogeological perspective, Yogyakarta is underlain by shallow aquifers with relatively high groundwater tables, as well as unconsolidated volcanic sand deposits derived from Mount Merapi eruptions. These conditions significantly increase susceptibility to pore water pressure buildup during strong ground shaking, thereby amplifying the likelihood of liquefaction.

In addition to strong ground shaking, earthquakes originating from the Opak Fault and other active fault systems in Central Java and Yogyakarta can trigger secondary geotechnical hazards such as liquefaction, especially in areas composed of unconsolidated alluvial and volcanic deposits (Mase et al., 2024; Mase et al., 2025). The 2006 Yogyakarta earthquake, for example, caused several liquefaction and lateral displacement events in the Bantul area and surrounding low-lying areas, indicating that the soil characteristics in this region are highly susceptible to increased pore pressure during seismic shaking. Recent studies have also emphasized that zones composed of young volcanic sediments and water-saturated sand layers common in the Yogyakarta Basin are highly susceptible to such failure mechanisms. Yogyakarta and Central Java, as well as other cities such as Bengkulu, Palu, which have been studied previously, exhibit high earthquake hazard potential due to active seismotectonic conditions, soft soil characteristics, and the presence of local faults such as the Opak Fault (Mase et al., 2024; Mase et al., 2025). The results of this study, which refine the spatial distribution of local seismicity and map active fault segments, are critical for identifying high-risk zones for liquefaction. The integration of these seismotectonic insights into liquefaction hazard modelling will improve the precision of the microzonation maps and provide a stronger basis for the design of earthquake-resistant structures, infrastructure planning, as well as community-based mitigation strategies aimed at reducing the cascading impacts of future seismic events in the Central Java and Yogyakarta regions. Collectively, these findings underscore that the interplay between subduction driven earthquakes, hydrogeological conditions, and local volcanic sediments plays a critical role in shaping liquefaction hazards. Therefore, integrated analyses combining seismotectonic data, groundwater conditions, and local geological characteristics are essential to improve the accuracy of hazard mapping and to strengthen disaster mitigation strategies in Indonesia’s subduction prone regions.

5. Conclusions

Significant activity is identified on the Baribis-Kendeng Fault, Kebumen-Meratus Trend (Central Java Fault), and Progo-Muria Lineament, each of which exhibits different fault mechanism characteristics and hypocenter depths varying from shallow, medium to deep. In addition, the distribution of seismicity around the Opak Fault Zone and associated structures (Oyo Fault, Ngalang Fault) confirms the existence of a complex fault network with potential inter-segment interaction and stress transfer, which has implications for multisegment earthquake scenarios. In the southern subduction zone of Java, the distribution of hypocenters shows a steep slab geometry, with a concentration of intermediate crustal earthquakes around 50 km deep and a cluster of intraslab earthquakes at 50–250 km depth. This activity is related to slab-pull forces and tectonic mechanisms in accretionary prisms and backthrust structures. The combination of nearshore shallow earthquakes, intermediate crustal earthquakes, and intraslab earthquakes broadens the hazard spectrum, including the potential for strong shaking, amplification of seismic waves due to thick volcanic deposits, local tsunamis in the event of seafloor deformation, and secondary hazards such as liquefaction in areas underlain by young alluvial deposits, volcanic sediments, and shallow aquifers. Moreover, several active faults that have not yet been officially mapped in the Indonesian Earthquake Source and Hazard Map (PUSGEN) may also pose a serious threat to the Central Java and Yogyakarta regions. Therefore, urgent attention and policy support from the government are required to update fault inventories, strengthen seismic hazard assessments, and integrate these findings into national disaster mitigation strategies. The findings confirm that seismic hazard modelling in Central Java and Yogyakarta should consider the multisource and multimechanism nature of earthquakes, with data integration into probabilistic models, detailed mapping of active earthquake sources, and simulation of shaking, liquefaction, and tsunami scenarios. This approach is needed to update risk maps, support hazard-based spatial planning, and improve infrastructure resilience and community preparedness.

References

Amin, T., N. Ratman, & S. Gafoer. (1999). Geological Map of Middle Part of Java [Map]. Geological Research and Department Centre.

Aribowo, S., Husson, L., Natawidjaja, D. H., Authemayou, C., Daryono, M. R., Puji, A. R., Valla, P. G., Pamumpuni, A., Wardhana, D. D., De Gelder, G., Djarwadi, D., & Lorcery, M. (2022). Active Back‐Arc Thrust in North West Java, Indonesia. Tectonics, 41(7), e2021TC007120. https://doi.org/10.1029/2021TC007120

Cilia, M. G., Mooney, W. D., & Nugroho, C. (2021). Field Insights and Analysis of the 2018 Mw 7.5 Palu, Indonesia Earthquake, Tsunami and Landslides. Pure and Applied Geophysics, 178(12), 4891–4920. https://doi.org/10.1007/s00024-021-02852-6

Damanik, R., Supendi, P., Widiyantoro, S., Rawlinson, N., Ardianto, A., Gunawan, E., Husni, Y. M., Zulfakriza, Z., Sahara, D. P., & Shiddiqi, H. A. (2021). Earthquake monitoring of the Baribis Fault near Jakarta, Indonesia, using borehole seismometers. Geoscience Letters, 8(1), 38. https://doi.org/10.1186/s40562-021-00209-4

Garini, S. A., Madlazim, & Rahmawati, E. (2014). Relokasi Hiposenter Gempa Bumi Di Sulawesi Tengah Dengan Menggunakan Metode Geiger Dan Coupled Velocity-Hypocenter. Jurnal Fisika, 3(2), 107–112. https://doi.org/10.26740/ifi.v3n2.p%25p

Gunawan, E., & Widiyantoro, S. (2019). Active tectonic deformation in Java, Indonesia inferred from a GPS-derived strain rate. Journal of Geodynamics, 123, 49–54. https://doi.org/10.1016/j.jog.2019.01.004

Hall, R., Clements, B., Smyth, H. R., & Cottam, M. A. (2007). A new interpretation of Java’s structure. https://archives.datapages.com/data/ipa_pdf/078/078001/pdfs/IPA07-G-035.pdf

Hayes, G. P., Wald, D. J., & Johnson, R. L. (2012). Slab1. 0: A three‐dimensional model of global subduction zone geometries. Journal of Geophysical Research: Solid Earth, 117(B1). https://doi.org/10.1029/2011JB008524

Hidayat, W., Sahara, D. P., Widiyantoro, S., Putu Raditya Ambara Putra, I., Shihab, N. H., Amalia, R., Harisandi, A., & Suharsono. (2021). Analysis of Hypocenter Relocation by Using Double Difference Method in a Deep Underground Mining. IOP Conference Series: Earth and Environmental Science, 873(1), 012099. https://doi.org/10.1088/1755-1315/873/1/012099

Hidayat, W., Sahara, D. P., Widiyantoro, S., Suharsono, S., Riyanto, E., Nukman, M., Wattimena, R. K., Melati, S., Sitorus, E., Nainggolan, T., & Putra, I. P. R. A. (2024). 4D time lapse tomography for monitoring cave propagation and stress distribution in Deep Mill Level Zone (DMLZ) PT Freeport Indonesia. Geomechanics and Geophysics for Geo-Energy and Geo-Resources, 10(1), 39. https://doi.org/10.1007/s40948-023-00718-w

Hidayat, W., Sahara, D. P., Widiyantoro, S., Suharsono, S., Wattimena, R. K., Melati, S., Putra, I. P. R. A., Prahastudhi, S., Sitorus, E., & Riyanto, E. (2022). Testing the Utilization of a Seismic Network Outside the Main Mining Facility Area for Expanding the Microseismic Monitoring Coverage in a Deep Block Caving. Applied Sciences, 12(14), 7265. https://doi.org/10.3390/app12147265

Hoffmann-Rothe, A., Ritter, O., & Haak, V. (2001). Magnetotelluric and geomagnetic modelling reveals zones of very high electrical conductivity in the upper crust of Central Java. Physics of the Earth and Planetary Interiors, 124(3–4), 131–151. https://doi.org/10.1016/S0031-9201(01)00196-0

Koulali, A., Mcclusky, S., Susilo, S., Leonard, Y., Cummins, P., Tregoning, P., Meilano, I., Efendi, J., & Wijanarto, A. B. (2017). The kinematics of crustal deformation in Java from GPS observations: Implications for fault slip partitioning. Earth and Planetary Science Letters, 458, 69–79. https://doi.org/10.1016/j.epsl.2016.10.039

Lallemand, S., Heuret, A., & Boutelier, D. (2005). On the relationships between slab dip, back‐arc stress, upper plate absolute motion, and crustal nature in subduction zones. Geochemistry, Geophysics, Geosystems, 6(9). https://doi.org/10.1029/2005GC000917

Liu, Y., Gu, Y., & Zhang, H. (2024). Seismic risk assessment and damage analysis: Emerging trends and new developments. Journal of Safety Science and Resilience, 5(3), 365–381. https://doi.org/10.1016/j.jnlssr.2024.04.005

Mase, L. Z., Amri, K., Ueda, K., Apriani, R., Utami, F., Tobita, T., & Likitlersuang, S. (2024). Geophysical investigation on the subsoil characteristics of the Dendam Tak Sudah Lake site in Bengkulu City, Indonesia. Acta Geophysica, 72(2), 893–913. https://doi.org/10.1007/s11600-023-01158-6

Mase, L. Z., Likitlersuang, S., Tobita, T., Chaiprakaikeow, S., & Soralump, S. (2020). Local site investigation of liquefied soils caused by earthquake in Northern Thailand. Journal of Earthquake Engineering, 24(7), 1181–1204.

Mase, L. Z., Refrizon, Pakpahan, N. S., Shelina, A., Wu, J., & Chaiyaput, S. (2024). Site Characteristics of Musi Fault Zone in Kepahiang Urban Area, Bengkulu Province, Indonesia. Transportation Infrastructure Geotechnology, 11(4), 2453–2482. https://doi.org//10.1007/s40515-024-00375-w

Mase, L. Z., Tanapalungkorn, W., Anussornrajkit, P., & Likitlersuang, S. (2025). Assessing liquefaction risk and hazard mapping in a high-seismic region: A case study of Bengkulu City, Indonesia. Natural Hazards, 121(6), 6597–6623. https://doi.org/10.1007/s11069-024-07057-3

Mase, L. Z., Tanapalungkorn, W., Ueda, K., & Likitlersuang, S. (2025). Non-linear Site Response Analysis of Liquefaction in the Port Area of Bengkulu City Due to Large Subduction Earthquakes. Transportation Infrastructure Geotechnology, 12(2), 87. https://doi.org/10.1007/s40515-025-00540-9

Muttaqy, F., Nugraha, A. D., Puspito, N. T., Sahara, D. P., Zulfakriza, Z., Rohadi, S., & Supendi, P. (2023). Double-difference earthquake relocation using waveform cross-correlation in Central and East Java, Indonesia. Geoscience Letters, 10(1), 5. https://doi.org/10.1186/s40562-022-00259-2

Novianto, A., Hidayat, W., & Bintarto, B. (2025). Interception of Meratus and Java Trends Structures as Controlling the Emergence of the Purwodadi Mud Volcano Complex. IOP Conference Series: Earth and Environmental Science, 1458(1), 012018. https://doi.org/10.1088/1755-1315/1458/1/012018

Novianto, A., Prasetyadi, C., & Hidayat, W. (2021). Mud Volcano: Revealing the Stratigraphy of Kendeng Basin, Indonesia. Open Journal of Yangtze Oil and Gas, 7(1), 48–64.

Pratama, C., Widjajanti, N., Pinasti, A., Ummi, R., Lestari, D., Sunantyo, T. A., Heliani, L. S., & Ulinnuha, H. (2021). A possible opak fault segment that caused the 2006 mw 6.3 yogyakarta earthquake and its future implication. 683(1), 012059.

Raharja, R., Ito, T., & Meilano, I. (2024). Evaluation of earthquake potential using a kinematic crustal block motion model in Java, Indonesia, based on GNSS observation. Journal of Asian Earth Sciences: X, 11, 100171. https://doi.org/10.1016/j.jaesx.2023.100171

Ramdhan, M., Priyobudi, P., Kristyawan, S., & Sembiring, A. S. (2020). Seismisitas di Wilayah Jawa Tengah dan Sekitarnya Berdasarkan Hasil Relokasi Hiposenter dari Empat Jaringan Seismik Menggunakan Model Kecepatan 3-D. Eksplorium, 41(1), 61. https://doi.org/10.17146/eksplorium.2020.41.1.5828

Saputra, H., Wahyudi, W., Suardi, I., Anggraini, A., & Suryanto, W. (2021). The waveform inversion of mainshock and aftershock data of the 2006 M6. 3 Yogyakarta earthquake. Geoscience Letters, 8(1), 9.

Satyana, A. H. (2008). Mud Diapirs and Mud Volcanoes in Depressions of Java to Madura: Origins, Natures, and Implications to Petroleum System. Proceedings of Indonesian Petroleum Association, Thirty-Second Annual Convention & Exhibition, 32, IPA08-G 139.

Soehaimi, A. (2008). Seismotektonik dan Potensi Kegempaan Wilayah Jawa. Indonesian Journal on Geoscience, 3(4), 227–240. https://doi.org/10.17014/ijog.vol3no4.20085

Susanta, F. F., Pratama, C., Aditya, T., Khomaini, A. F., & Abdillah, H. W. K. (2019). Geovisual Analytics of Spatio-Temporal Earthquake Data in Indonesia. JGISE: Journal of Geospatial Information Science and Engineering, 2(2). https://doi.org/10.22146/jgise.51131

Syahbana, A. J., Iqbal, P., Irsyam, M., Asrurifak, M., & Hendriyawan, H. (2021). Smoothed gridded seismicity effect for land-use development, case study: Kalimantan Island, Indonesia. Rudarsko-Geološko-Naftni Zbornik, 36(3).

Tsuji, T., Yamamoto, K., Matsuoka, T., Yamada, Y., Onishi, K., Bahar, A., Meilano, I., & Abidin, H. Z. (2009). Earthquake fault of the 26 May 2006 Yogyakarta earthquake observed by SAR interferometry. Earth, Planets and Space, 61(7), e29–e32. https://doi.org/10.1186/BF03353189

Um, J., & Thurber, C. (1987). A Fast Algorithm for Two-Point Seismic Ray Tracing. Bulletin of the Seismological Society of America, 77(3), 972–986. https://doi.org/10.1785/BSSA0770030972

Upomo, T. C., Chang, M., Kusumawardani, R., Prayitno, G. A., Kuo, C.-P., & Nugroho, U. (2023). Assessment of Petobo Flowslide Induced by Soil Liquefaction during 2018 Palu–Donggala Indonesian Earthquake. Sustainability, 15(6), 5371. https://doi.org/10.3390/su15065371

Wagner, D., Koulakov, I., Rabbel, W., Luehr, B.-G., Wittwer, A., Kopp, H., Bohm, M., Asch, G., & MERAMEX scientists. (2007). Joint inversion of active and passive seismic data in Central Java. Geophysical Journal International, 170(2), 923–932.

Waldhauser, F. (2001). hypoDD -- A Program to Compute Double-Difference Hypocenter Locations (hypoDD version 1.0—03/2001).

Waldhauser, F., & Ellsworth, W. L. (2000). A double-difference earthquake location algorithm: Method and application to the northern Hayward fault, California. Bulletin of the Seismological Society of America, 90(6), 1353–1368.

Xu, S., Fukuyama, E., Ben-Zion, Y., & Ampuero, J.-P. (2015). Dynamic rupture activation of backthrust fault branching. Tectonophysics, 644, 161–183.

Zhang, H. (2003). Double-Difference Tomography: The Method and Its Application to the Hayward Fault, California. Bulletin of the Seismological Society of America, 93(5), 1875–1889. https://doi.org/10.1785/0120020190

Zhang, H., & Thurber, C. (2006). Development and Applications of Double-difference Seismic Tomography. Pure and Applied Geophysics, 163(2–3), 373–403. https://doi.org/10.1007/s00024-005-0021-y

Zulfakriza, Z., Nugraha, A. D., Heryandoko, N., Ry, R. V., Muttaqy, F., Andika, A., Azhari, M. F., Putra, A. S., Palgunadi, K. H., Cummins, P. R., Supendi, P., Lesmana, A., Sahara, D. P., & Puspito, N. T. (2024). Seismic source analysis of the destructive earthquake November 21, 2022, Mw 5.6 Cianjur (Indonesia) from relocated aftershock. Scientific Reports, 14(1), 12142. https://doi.org/10.1038/s41598-024-60408-9

Funding

This research was supported by funding from the Institute for Research and Community Service (LPPM), UPN “Veteran” Yogyakarta, under Grant No. 422/UN62.21/PG.00.00/2025, awarded to Wahyu Hidayat.

Author’s contribution

Wahyu Hidayat: conceptualization, supervision, methodology, writing and interpretation, Nanda Agustina validation, processing and visualitation, Ekalia Ardenti data curation, software, resources supervision, validation, Suharsono original draft and writing, review and editing. Ardian Novianto project administration and review and editing. Uli Ulfa review and editing. Maman Rohaman software, resources supervision, Herry Riswandi writing and interpretation.

All authors have read and agreed to the published version of the manuscript.