STRESS TRANSFER AND SEISMICITY ALONG THE SOUTHERN GREAT SUMATRAN FAULT: INSIGHTS FROM THE 1994 LIWA EARTHQUAKE

1. Introduction

On February 15, 1994, at 17:08 UTC, a Mw 6.8 earthquake struck west of Lampung in southern Sumatra, Indonesia (see the yellow star in Figure 1). This destructive mainshock occurred on the southern segment of the Great Sumatran Fault (GSF) near the city of Liwa, in the Semangko Fault. Around 207 people died, over 2,000 were injured, and 75,000 were displaced; most of the damage and fatalities occurred in the Liwa area. The earthquake damaged over 6,000 homes, shops, and government buildings (Widiwijayanti et al., 1996). Since the early 20th century, at least three significant earthquakes (M ≥ 6.0) have occurred around Liwa (red star in Figure 1): (1) the 1908 event on the western Semangko segment of the GSF; (2) the 1933 event on the Kumering segment, Ms 7.5 (Natawidjaja & Triyoso, 2007); and (3) the 1994 Liwa earthquake on the Semangko segment, Mw 6.8 (Soehaimi et al., 2002).

According to Bellier et al. (1997), Natawidjaja et al. (2017), and Wang et al. (2023), the GSF is about 1,900 km long and an important part of the strike-slip fault system that allows the Indo-Australian–Eurasian plate to move at an angle. Regional convergence rates are ~44 mm/yr in the Andamans, ~60 mm/yr in southern Sumatra, and ~68 mm/yr in southern Java (DeMets et al., 2010). Based on these rates, the GSF can handle several millimeters of trench-parallel motion per year. Interseismic geodetic modelling puts the system-wide slip at about 15 mm/yr. The slip rates in the southern GSF are about 18 mm/yr on the Manna segment and about 12 mm/yr on Kumering. Musi is less well known, with locking depths of about 1–20 km (Lubis et al., 2024). The Semangko section divides into the West and East Semangko Faults close to Liwa, moving at rates of 16.5 ± 2 and 12.5 ± 2 mm/yr, respectively (Alif et al., 2020). Even though these rates are significant, no Mw ≥ 6.0 earthquake has occurred on the southern GSF since the 1994 Liwa event, indicating continuous stress accumulation and unresolved seismic potential (Salman et al., 2020; Hutchings & Mooney, 2021; Triyoso & Suwondo, 2023; Kurniawan et al., 2024).

Coulomb Failure Stress (ΔCFS) modelling has been widely applied to quantify stress changes caused by earthquakes and ongoing loading, and is effective for identifying regions with increased failure potential. Numerous studies have demonstrated that static stress changes on the order of 0.01–0.1 MPa can promote or inhibit fault slip, depending on fault orientation and pre-existing conditions (King & Devès, 2015; Pope & Mooney, 2020; Segou & Parsons, 2020; Mildon et al., 2022). ΔCFS modelling is now used in many different types of tectonic settings, including strike-slip systems in California (Pope & Mooney, 2020), thrust faulting in the Himalayas, India (Parija et al., 2021), and normal-fault systems in the Central Apennines, Italy (Valentini et al., 2024). In Indonesia, ΔCFS evaluations have addressed segments of the GSF following megathrust earthquakes along the Sunda subduction zone (Rafie et al., 2023; Sukrungsri et al., 2024) and several onshore sequences, such as the Pasaman, West Sumatra sequence (Supendi et al., 2023; Wulandari et al., 2023) and the Angkola segment, North Sumatra (Sahara et al., 2018). However, applications of ΔCFS are rare in southern Sumatra. To close this gap, we need to do a thorough job of characterizing ΔCFS patterns. This will help us better assess seismic hazards, find areas where seismic energy is more likely to be released, prioritize mitigation efforts, and set requirements for infrastructure resilience.

Triyoso and Suwondo (2023) recently found high compression around Liwa along the Kumering–Semangko section of the GSF. They also found high slip rates and a high correlation dimension, all of which are signs of high stress. Since Liwa is close to the GSF and could make it stronger, we need an integrated risk framework that clearly limits the crustal GSF source along with megathrust and intermediate-depth sources. We use the double-difference method to relocate earthquakes more accurately and combine these results with static ΔCFS modelling to (i) figure out how much stress the 1994 mainshock changed things; (ii) see how well the relocated clusters match up with the ΔCFS lobes; and (iii) see what this means for future rupture potential and seismic risk along the southern GSF. This analysis will help discriminate whether the 1994 rupture involved a previously unmapped GSF segment or a distinct, as-yet unidentified fault, thereby tightening physics-based constraints for regional hazard models.

Southern Sumatra remains under-documented despite pronounced tectonic complexity; pre-2010 phase picks are not publicly available and station coverage is heterogeneous, which constrains relocation fidelity before 2010. Our hypoDD workflow therefore uses 1-D velocity models and differential times to reduce model dependence, and the ΔCFS analysis adopts static Okada solutions in a homogeneous, isotropic half-space (Okada, 1992). These assumptions provide a first-order, physics-based baseline for the southern GSF but do not capture viscoelastic relaxation or afterslip. This gap motivates targeted acquisition (denser seismic arrays), integration of GNSS/InSAR and 3-D tomography, and focused paleoseismic work on ESF-A/WSF-B to refine fault geometry, locking depth, and recurrence in future studies.

Figure 1. (a) Sumatra in the Sunda subduction system. The red triangles represent the Sunda Trench, and the arrows show the trench-normal convergence rates, which are about 53–61 mm yr⁻¹. The Barisan Mountains extend from Banda Aceh to the Sunda Strait via the Great Sumatran Fault (GSF; red lines). The blue circles show where most people live, and the green triangles show where BMKG stations are located. The study area is shown by the black rectangle, and the inner rectangle shows the borders of the Indo-Australian, Eurasian, and Pacific plates. (b) Zoom in on the area of study. The GSF separates into the Kumering Fault and the West and East Semangko Faults (WSF and ESF). The yellow star shows the Mw 6.8 Liwa mainshock on February 15, 1994, with a focal mechanism. The red stars show the 1933 Ms 7.5 Kumering and 1908 occurrences on WSF. The main stress frame is shown with SHmax ∼ N30°E and Shmin ∼ N294°E, aligning with initial motion and regional stress parameters (e.g. Widiwijayanti et al., 1996; Heidbach et al., 2018).

2. Seismotectonics of the Great Sumatran Fault (GSF)

Sumatra Island is located on the southern edge of the Sunda Plate, which engages morphologically with the subducting Indo-Australian Plate beneath it. The convergence component of about 50–70 mm/year has been divided up: the part that runs parallel to the trench is mostly taken up by the right-lateral strike-slip fault system of the GSF, and the part that runs perpendicular to the trench is taken up by the Sunda megathrust (Sieh & Natawidjaja, 2000). The GSF extends from Banda Aceh in the north, across the Barisan Mountains, and ends at the Sunda Strait (Figure 1a; Natawidjaja, 2018). It is one of the most active fault systems in Indonesia. The GSF was first divided into 19 segments with various characteristics (Sieh & Natawidjaja, 2000). Later, PUSGEN (Pusat Studi Gempa Nasional) improved the map by adding about 40 additional segments (Irsyam et al., 2017).

A series of large megathrust earthquakes off the coast of Sumatra have raised the seismic risk along the GSF. These earthquakes have transferred static stress and are likely to put more strain on the GSF (Rafie et al., 2023). The Mw 9.2 Sumatra-Andaman earthquake on December 26, 2004, the Mw 8.6 Nias-Simeulue earthquake on March 28, 2005, and the Mw 8.4 Bengkulu earthquake on September 12, 2007, all came first in this sequence. Because of these major events, stress accumulation along the GSF is believed to have increased, thereby elevating the likelihood of subsequent earthquakes along this fault system. Rafie et al. (2023) emphasize the importance of continuous seismic monitoring and preparedness in the region, as static stress transfer can potentially trigger new seismic activity with significant impacts on communities and infrastructure throughout Sumatra.

This effect is reflected in the growing frequency of moderate (M~6) earthquakes along the GSF after these large megathrust events. Consistently, a series of moderate earthquakes have repeatedly occurred along the fault, including two Mw 6.4 and 6.3 events near Lake Singkarak in 2007 (Daryono et al., 2012), the Mw 6.0 Toru earthquake in 2008 that ruptured the Toru segment, the Mw 6.6 Dikit seismic gap earthquake in 2009 (Hurukawa et al., 2014), the Mw 6.1 creeping segment earthquake in Aceh in January 2013 (Ito et al., 2016), and the Mw 6.1 Pasaman earthquake in February 2022 involving the Angkola, Sumpur, Sianok, and Barumun segments (Supendi et al., 2023; Wulandari et al., 2023). This sequence underscores the dual influence of both the megathrust and the GSF fault systems in shaping seismic hazard patterns across Sumatra.

In the southern sector, the Semangko segment forms a branching fault system, the West Semangko Fault (WSF), East Semangko Fault A (ESF-A) and East Semangko Fault B (ESF-B), which are separated by a step-over structure in the Ranau Caldera region. This geometry causes the deformation regime to transition from predominantly strike-slip to transtensional within the pull-apart basin. The presence of Quaternary volcanic-sedimentary infill overlying the Barisan Mountains arc basement produces a stiffness contrast and weak zones, focusing strain along the basin margins and minor splays. Regionally, the obliquity of the Sunda subduction arranges a transpressional regime characterized by right-lateral strike-slip accumulation along the GSF. However, at the Ranau step-over, stress field rotation introduces an oblique-normal faulting component on secondary faults, consistent with shallow crustal seismicity clusters near Liwa (Heidbach et al., 2018; Triyoso & Suwondo, 2023). P-wave first-motion mechanism inversion reveals a strike-slip regime with a maximum horizontal stress axis of approximately N30°E and a minimum of about N294°E (e.g. Widiwijayanti et al., 1996; Heidbach et al., 2018). Concurrently, heterogeneity in slip rates between branches and variations in locking depth generate a differential stress gradient that increases the likelihood of cross-segment triggering (Alif et al., 2020; Lubis et al., 2024). Post seismic static stress transfer following the 1994 Liwa earthquake (ΔCFS) is also thought to have reinforced local loading on splays and the caldera margins. The combination of step-over fault geometry, weaker basin materials, and interconnected fault sources makes the southern Semangko–Ranau zone a critical focus for seismic hazard assessment.

Table 1. The source parameters of the Liwa earthquake by Widiwijayanti et al. (1996) and Soehaimi et al. (2002). The length and the width were estimated based on the empirical relations of Wells & Coppersmith (1994).

3. Data and Methods

The earthquake dataset covers a span of 30 years from February 1994 to January 2024 (following the 1994 Mw 6.8 Liwa earthquake) and is sourced from BMKG's permanent seismic network in Lampung. Within this timeframe, we compiled 498 events (Mw 0.9–5.0) monitored by 209 BMKG stations in the southern Sumatra and western Java regions (see Figure 1). For the period 1994–2009, only catalog data are available because of limitations in station density and phase archives. From 2010 to 2024, P- and S-wave arrival times are available, enabling double-difference relocation. Operationally, the initial hypocenters were determined through linear inversion in SeisComP3 (Supendi et al., 2022). This initial solution is sufficient for monitoring, but accuracy may suffer due to the subsurface velocity model's simplification and the uneven distribution of propagation paths, particularly in complex tectonic regions. Active faults around the Liwa earthquake source were compiled from PUSGEN (Irsyam et al., 2017). The direction of the fault plane and the type of displacement (focal mechanism) of the 1994 Mw 6.8 main earthquake are based on the research based on Soehaimi et al. (2002). The aftershock distribution was taken from the seismic observations of Widiwijayanti et al. (1996) and digitized to get the event coordinates, limit the shape of the rupture, and make sure the static stress change modelling was correct. The source dimensions and slip magnitude were estimated using the well-established empirical relationship developed by Wells & Coppersmith (1994), known for its accuracy in seismic source parameter estimation, and are detailed in Table 1.

We relocated hypocenters using the hypoDD double-difference algorithm (Waldhauser & Ellsworth, 2000; Waldhauser, 2001), which reduces the difference in travel time between two nearby events as much as possible. The main idea is that when the distance between events is much smaller than the distance between sources and stations, the paths of the rays are almost the same, and they sample similar velocity structures. This means that differential times can be used to get very accurate relative locations. We did two separate inversions with the same parameters to see how the 1-D velocity model affected the accuracy of the relocation (see Figure 2): CRUST1.0 (Laske et al., 2013) and AK135 (Kennett et al., 1995; Hu et al., 2023; Pranata et al., 2023). The quality of the move was judged by looking at the distribution of differential-time residuals, which were shown as histograms with the same number of bins and axes for both cases. The reference model was selected based on a stronger concentration of residuals near 0 s (a sharper peak), indicating a better overall data–model fit.

We limited the distance between events to 100 km, let up to 30 nearest neighbors per event, and needed at least 10 differential-time links to ensure the network was connected. To keep differential-time quality, we also limited the distance between the station and the cluster centroid to 300 km. These constraints were applied uniformly to both velocity-model tests to enable a fair comparison. We chose the parameters by running repeated sensitivity tests until we found a stable configuration that moved things consistently across a range of possible inputs. This workflow mitigates bias from unmodelled lateral velocity heterogeneity and enhances the spatial coherence of seismic clusters, thereby improving the interpretability of subsequent static stress modelling in structurally complex fault systems (e.g. Muttaqy et al., 2023; Supendi et al., 2023).

Figure 2. One-dimensional (1D) velocity models tested in the relocation experiments. Solid lines show the CRUST1.0 model and dashed lines the AK135 model; blue curves denote P-wave velocity (Vp) and red curves S-wave velocity (Vs). Depth is in kilometers, positive downward.

We looked at how stress changed after the earthquake using static ΔCFS modelling with Coulomb 3.3 (Toda et al., 2011), which was based on the work of King et al. (1994). We define

ΔCFS = Δτ - μ′ Δσ_n (1)

Where Δτ is the change in shear stress, Δσₙ is the change in normal stress (positive for unclamping), and μ′ is the effective coefficient of friction. We calculated stress changes on two types of receiver faults: (i) planes that were perfectly aligned with the Semangko Fault's local strike at depths of 5, 10, 15, and 20 km (these were the seismogenic layers); and (ii) mapped active faults around the mainshock that were used as receiver faults and had their own strike–dip–rake. These faults came from the national active-fault database (Irsyam et al., 2017). Receiver faults are potential slip surfaces (either hypothetical, optimally oriented planes or mapped structures) on which we evaluate whether stress changes would promote (ΔCFS > 0) or inhibit (ΔCFS < 0) failure. The mainshock was shown as a flat displacement whose size, shape, rake, strike, and dip were limited by the focal-mechanism solution and the geometry of the aftershock. Wells and Coppersmith's (1994) empirical relations were used to estimate the fault parameters. We assumed a homogeneous, isotropic elastic half-space with Poisson’s ratio = 0.25 (Xu et al., 2016). To bracket frictional uncertainty, simulations were run for μ′ = 0.2, 0.4, 0.6, and 0.8 (Wulandari et al., 2023). To make sure the numbers were correct, CFS values were taken at the hypocenters and at random locations inside the convex hull of events that were widened by 0.2°. Interpolation was linear (nearest neighbor checked for robustness). We report the fractions of ΔCFS > 0 and ΔCFS > 0.01 MPa, median values, and results of a one-sided Kolmogorov–Smirnov test comparing the two populations.

We measured the correlation between ΔCFS and relocated seismicity, in addition to spatial coincidence. For each hypocenter, ΔCFS was collected using depth-specific grids (5-20 km) estimated for the 1994 Liwa source and an effective friction range (μ′ = 0.2-0.8). For a fair baseline, N=5000 random sites were chosen inside the convex hull of the event distribution extended by a 0.2° buffer, and ΔCFS was derived by linear interpolation (nearest checked for robustness). We quantified the association using the fraction of events with ΔCFS > 0 and ΔCFS ≥ 0.01 MPa (10 kPa), median contrasts, and a one-sided Kolmogorov–Smirnov test (H₁: events exhibit higher ΔCFS than random points). Fraction-above-threshold metrics are commonly used because the random domain is largely dominated by ΔCFS = 0 outside stress-change zones (e.g. King et al., 1994; Toda et al., 2005), thereby complementing whole-distribution tests. For the Liwa mainshock, Coulomb stress modelling was performed under the assumption of a static, elastic, homogeneous half-space based on the Okada dislocation formulation (Okada, 1992), with the rupture represented as a planar fault geometry. The receiver faults were defined by mapped strands of the ESF and WSF. Sensitivity analyses were conducted to evaluate the influence of effective friction (μ′), source depth, interpolation scheme, and the extent of the random-domain buffer.

4. Results

4.1. Seismicity Relocation Quality

Using identical inversion parameters, we evaluated 16,000 combined P- and S-wave arrivals. The post-relocation residual distributions show a clear contrast between the two 1-D velocity models (see Figure 3). CRUST1.0 yields a tight, nearly symmetric distribution centered at 0 s, with most residuals within ± 0.3–0.4 s. By contrast, AK135 produces a broader, heavy-tailed distribution extending to ± 0.7–0.8 s and a slight negative bias, implying an underestimation of crustal velocities along local paths. The sharper near-zero peak and reduced spread achieved with CRUST1.0 indicate a superior data–model fit and a more realistic representation of shallow crustal heterogeneity in the study area. We therefore adopted CRUST1.0 for the final relocated catalog and all subsequent analyses. Since station coverage and phase archives are limited before 2010, events from 1994–2009 are treated as background, whereas 2010–2024 earthquakes with phase picks were relocated. Figure 4 contrasts the pre-relocation catalogue with the double-difference solution. The operational locations exhibit diffuse epicenters and a pronounced fixed-depth band near at ~10 km in the A–B section, typical of 1-D linearized inversions with uneven path coverage. The hypoDD relocation markedly reduces depth scatter and collapses epicenters into coherent clusters that delineate the Semangko Fault system, with prominent concentrations beneath Liwa and Semangko Bay. The improved geometry provides a more reliable basis for subsequent ΔCFS analysis.

Figure 3. Post-relocation differential travel-time residuals from hypoDD for two 1-D velocity models. Histograms use identical axes and binning. CRUST1.0 (green) exhibits sharp, near-zero peaks with a narrow spread (≈±0.3–0.4 s), whereas AK135 (orange) is broader with heavier tails (±0.7–0.8 s) and a slight negative bias. The tighter, near-zero-centered distribution indicates a superior data–model fit; accordingly, CRUST1.0 is adopted for the final relocations and subsequent analyses.

Figure 4. (a) Operational catalogue locations for events used in this study; (b) relocations obtained with hypoDD. Symbol size scales with magnitude; colour indicates hypocentral elevation (km). The dashed A–B line marks the profile trace; lower panels show A–B cross-sections with surface topography (Tozer et al., 2019; Wessel et al., 2019). Of the 498 events with usable picks, 416 were retained after relocation; removed events failed the graph-connectivity/quality criteria (≥10 differential links, ≤100 km inter-event distance, ≤300 km station-cluster offset). Relocation suppresses the artificial fixed-depth band near ~10 km and tightens spatial clustering near Liwa, Semangko Bay, and along the Semangko Fault system

4.3. Statistical Validation of ΔCFS and Seismicity

To examine the relationship between stress perturbations and relocated seismicity, we performed a straightforward statistical assessment. Each relocated event from 2010 to 2024 was assigned a ΔCFS value at its hypocentral location using the static Coulomb stress-change model. The dataset was then divided into two categories: positive ΔCFS (>0.0 MPa) and negative ΔCFS (<0.0 MPa). The results indicate that approximately 72% of relocated earthquakes (n = 345 of 480) occurred within regions of positive ΔCFS, whereas about 28% (n ≈ 135) were in negative stress-change domains.

This asymmetric distribution supports the hypothesis that static stress redistribution from the 1994 Liwa earthquake promoted subsequent seismicity on adjacent fault segments. A Kolmogorov–Smirnov (K–S) test further confirms that the observed seismicity distribution is statistically distinct (p < 0.05) from a uniform random distribution with respect to ΔCFS. This strengthens the interpretation of the positive lobes of the Coulomb stress transfer function as preferential sites for earthquake triggering.

According to ECDF and histogram analyses, relocated events exhibit higher ΔCFS values than randomly sampled sites within the convex-hull-plus-buffer domain. Seventy-four percent of the events occurred at locations with ΔCFS > 0, while 73% exceeded the commonly applied promotion threshold of 10 kPa (0.01 MPa). The median ΔCFS at event sites (0.108 MPa) is substantially greater than the random-site median (0.002 MPa). A one-sided Kolmogorov–Smirnov test yielded D = 0.051 and p = 0.128, indicating a rightward shift in the event distribution but reduced statistical power due to the clustering of random samples near ΔCFS ≈ 0. Taken together, these observations highlight a pronounced preference for positively loaded lobes in post-1994 seismicity.

Figure 5. Spatial distribution of static ΔCFS changes at 10 km depth following the 1994 Mw 6.8 Liwa earthquake. Panels (a)–(d) show ΔCFS for μ′ = 0.2, 0.4, 0.6, and 0.8, respectively. Red and blue lobes indicate regions of positive and negative ΔCFS, respectively, with values ranging from –0.10 to +0.10 MPa. Aftershock locations (yellow–green circles) and mapped faults are overlaid for reference.

5. Discussion

5.1. Earthquake Relocation

Before relocation, the hypocenter distribution appeared more scattered and unstructured, with numerous events concentrated at fixed depths. This pattern reflects the limitations of conventional initial location methods. After relocation (right panel), the hypocenter distribution became significantly more focused, forming geologically plausible clusters, particularly around Semangko Bay and the WSF and ESF-A segments. These results indicate that the double difference method effectively mitigates the influence of local velocity heterogeneities and systematic inter-station errors, consistent with findings by Waldhauser (2001) and Supendi et al. (2023). Post-relocation results reveal a pronounced linear alignment of seismic events, highlighting a clear spatial coherence between the relocated hypocenters and known surface fault traces. Notably, the WSF and ESF-A segments exhibit intensified clustering of seismicity that follows the orientation of the mapped fault lines, reinforcing their seismotectonic significance. In addition, a substantial increase in seismic event density is observed within the central portion of the study area. This heightened concentration may indicate zones of elevated stress accumulation or the activity of previously unrecognized secondary structures. The spatial distribution of these events suggests the potential involvement of complex fault geometries, including antithetic fault systems, graben structures, or fault bifurcations, in accommodating crustal deformation. These observations are consistent with previous studies in tectonically active volcanic regions, where intricate subsurface fault networks often govern the pattern and intensity of seismicity (Supendi et al., 2022). The improved resolution achieved through the double difference relocation method thus offers valuable insight into the underlying fault architecture and may aid in the identification of seismogenic structures that remain undetected through surface mapping alone.

The depth cross-section reveals a significant contrast in hypocenter distribution before and after relocation. Prior to relocation, many events were artificially constrained at a uniform depth near 10 km, which likely reflects the limitations of conventional location algorithms rather than the actual seismogenic structure. Such vertical clustering can obscure interpretations of fault geometry and crustal deformation processes. After applying the double difference relocation, the events display a more geologically plausible vertical spread, ranging from 5 to 20 km in depth. This broader distribution delineates an active seismogenic zone and agrees with regional tomographic results and crustal velocity models reported by Hu et al. (2023). Most of the seismicity concentrates between 8 and 15 km depth, coinciding with the transitional zone from brittle to ductile deformation regimes. This observation supports the hypothesis that mid-crustal earthquakes in the GSF region are associated with thermal and rheological gradients, consistent with the findings of Pranata et al. (2023). Overall, the relocated seismicity provides improved constraints on the depth extent of faulting and enhances our understanding of the mechanical layering within the crust.

Several earthquake clusters are observed outside the main fault traces, particularly in the southwestern sector of the study area. These clusters are not spatially aligned with any documented surface faults, suggesting the possible existence of buttress structures, blind faults, or antithetic faulting. Similar patterns have been noted by Supendi et al. (2023), who associated certain clusters in the Semangko region with unmapped minor fault segments not reflected in current regional geological maps.

Figure 6. Spatial distribution of Coulomb Failure Stress (ΔCFS) changes induced by the 1994 Mw 6.8 Liwa earthquake, calculated at 10 km depth with an effective friction coefficient (μ′) of 0.4. The four-lobed stress pattern is characteristic of strike-slip fault rupture, showing areas of stress increase (red) and decrease (blue). Relocated seismicity from 2010–2024 is overlaid, with symbol sizes indicating event magnitude and colours representing focal depths. Three primary seismic clusters (A, B, and C) are delineated by dashed yellow ellipses.

5.2. Coulomb Failure Stress Distribution and Seismotectonic Implications

The ΔCFS distribution model presented in this study exhibits a consistent pattern of positive stress changes along segments of the Semangko Fault near the source area of the 1994 Liwa earthquake. The presence of elongated positive stress lobes (in red), aligned parallel to the fault and asymmetrically distributed on both sides of the rupture, indicates the occurrence of static stress transfer to adjacent fault segments with spatial and structural proximity (see Figure 6). Theoretically, this pattern supports the potential for future triggered seismicity, as originally proposed by King et al. (1994) and further substantiated by subsequent studies (Toda et al., 2005; Pope and Mooney, 2020).

Validation of the ΔCFS model against the historical seismicity distribution reveals a strong spatial correlation, particularly within clusters A, B, and C, which exhibit increased post-1994 earthquake activity (see Figure 6). These clusters coincide with areas of modelled stress increase exceeding 0.05 bar, supporting the interpretation that static stress redistribution enhanced the likelihood of subsequent seismicity in these zones. Similar findings have been reported by Xiong et al. (2017) and Mildon et al. (2018), who demonstrated that seismicity concentrations often align with positive ΔCFS anomalies in transfer zones and fault segment intersections.

The ΔCFS model is consistent with the findings of Supendi et al. (2022), who identified the presence of a graben system and antithetic faults beneath Semangko Bay features that contribute to the complex seismic response observed in the area. Activation of antithetic faults may locally amplify ΔCFS through mechanical interactions and the internal structural geometry of the graben. This highlights the importance of incorporating minor fault elements in static stress transfer modelling to adequately capture the complexity of stress redistribution. At the regional scale, these results also align with the findings of Sahara et al. (2018), which demonstrated that stress heterogeneity within bifurcation zones of the GSF such as the Barumun, Angkola, Toru segment can lead to the concentration of moderate-to-large magnitude earthquakes. Bifurcation and step-over zones along the GSF, as previously characterized in paleoseismic studies (Bellier et al., 1997), suggest that fault geometry complexity significantly enhances local stress accumulation and contributes to the occurrence of seismic clusters. In this case, fraction-based diagnostics (ΔCFS > 0 and ΔCFS ≥ 10 kPa) provide stronger evidence for static triggering than whole-distribution tests, with ~74% of events located in positively stressed regions. This proportion is consistent with previously reported ranges (≈60–80%) for other strike-slip fault systems.

Figure 7. Three-dimensional visualization of ΔCFS distribution for right-lateral slip across segmented faults in the southern Sumatra region. Blue is the rupture plane of the Liwa earthquake. The model simulates stress changes following the 1994 Liwa earthquake, highlighting key fault segments including the Kumering Fault, West and East Semangko Faults (A and B), Ujung Kulon Fault, and Semangko Graben.

Figure 8. Coulomb failure stress changes induced by the 1994 Mw 6.8 Liwa earthquake, modelled at four different depths: (a) 5 km, (b) 10 km, (c) 15 km, and (d) 20 km. The stress lobes show systematic propagation along and across the mapped fault segments, with notable positive ΔCFS regions coinciding with the WSF-B, ESF-A, and Kumering Fault.

5.3. Three-Dimensional Stress Transfer and Fault Interaction in Southern Sumatra

The model indicates that several fault segments, including ESF-A, WSF, and the Kumering Fault, experienced an increase in ΔCFS following the 1994 mainshock (see Figure 7). Areas of positive ΔCFS are represented by red shading on the grid surface, theoretically corresponding to zones brought closer to failure and thus exhibiting an elevated potential for triggered seismicity. In contrast, regions of negative ΔCFS (blue) are predominantly observed along the rupture segment of the 1994 Liwa earthquake, suggesting a stress decrease in those areas.

The elevated ΔCFS observed on the ESF-A and ESF-B segments corresponds with their fault geometry being aligned or subparallel to the rupture plane of the 1994 mainshock (see Figure 7). This observation is consistent with findings by Parija et al. (2021) and Supendi et al. (2022), which suggest that fault segments with orientations optimally aligned with the regional stress field tend to receive greater static stress transfer, thereby increasing their susceptibility to secondary seismic energy release.

Furthermore, the existence of minor structures and fault geometry complexity, such as the East Semangko Fault A/B, Ujung Kulon Fault, and Semangko Graben, emphasizes the significant role of local structural features in static stress redistribution. The 3D ΔCFS model (see Figure 7) supports Bellier et al. (2006)'s findings and Sahara et al. (2018)'s en-echelon fault patterns inside the GSF bifurcation zone, revealing positive stress concentrations that propagate southward and eastward by minor faults. These results show that minor faults cannot be ignored in tectonically complex systems such as southern Sumatra. Integrating ΔCFS modelling with high-resolution hypocenter relocation improves the understanding of stress evolution in fault networks, as seen in other complex locations including California, the Himalayas and Taiwan (Pope & Mooney, 2020; Parija et al., 2021; Ishimaru et al., 2025).

The spatial distribution of ΔCFS across depths of 5–20 km (see Figure 8) reveals patterns that reflect mechanical heterogeneity within the crust. At shallower depths (5–10 km), stress changes are more laterally dispersed, whereas at intermediate depths (15–20 km), ΔCFS becomes more concentrated along the trace of the main fault. The plots indicate that relocated events preferentially occur in ΔCFS-positive regions. This depth-dependent variation is most likely caused by variations in effective friction and rheological characteristics with increasing depth (Cattin et al., 2009), as well as the distribution from previous slip events, for instance, the 1933 and 1994 earthquakes.

The ΔCFS results presented here exhibit features comparable to global case studies, such as those in the California region (Pope & Mooney, 2020) and the Western Himalayas (Parija et al., 2021), particularly in the asymmetric distribution of stress lobes and the potential for post-seismic seismicity clustering. However, Sumatra’s tectonic setting is distinct, primarily governed by forearc sliver motion, which plays a critical role in accelerating stress accumulation along the GSF (Natawidjaja & Triyoso, 2007; Rafie et al., 2021). This tectonic control is further emphasized by Hutchings and Mooney (2021), who reported that the rate of stress accumulation in Sumatra is among the highest observed in subduction zones worldwide.

Figure 9. Joint diagnostics of ΔCFS at relocated earthquake hypocenters (orange) versus random samples (blue) drawn from within the convex-hull-plus-buffer domain. (Top) Empirical cumulative distribution functions (ECDFs) with vertical reference lines at ΔCFS = 0 and ΔCFS = 0.01 MPa (10 kPa). (Bottom) Histograms constructed using a Freedman–Diaconis bin width. Annotated statistics include Frac (ΔCFS > 0) = 0.74, Frac (ΔCFS ≥ 0.01 MPa) = 0.73, median(event) = 0.108 MPa, median(random) = 0.002 MPa, KS D = 0.051, and p = 0.128. Collectively, the plots indicate that relocated events preferentially occur in ΔCFS-positive regions relative to random expectation.

The model results, particularly for the Semangko segment, support concerns regarding the accumulation of unreleased stress since the last major earthquake in 1994. The western section of this segment displays positive ΔCFS values exceeding 0.1 bar across all modelled depths. According to the time-evolution ΔCFS model by Rafie et al. (2023), a fault segment with a stress accumulation rate of 17 kPa/year would require approximately 100–150 years to reach a critical seismic threshold. This suggests that the Semangko and Kumering segments would be ideal candidates for a future large-magnitude event if stress is not released by seismic creep or small earthquakes.

The presence of elevated ΔCFS values in this model further indicates that the Semangko Fault system is currently in a critical stress state, particularly in the southwestern section, which has not experienced significant stress release since the 1933 earthquake. This observation is broadly consistent with the concept of “stress recovery time” along the GSF, which is estimated to range between 100 and 200 years for stress levels to return to failure thresholds following a major rupture. These findings carry important implications for seismic hazard modelling in southern Sumatra. In particular, the potential for stress concentration within transition zones between fault segments (e.g. between the Semangko and Ujung Kulon faults) should be considered in future hazard zoning updates (Irsyam et al., 2017). Moreover, this study highlights the need for continuous monitoring of microseismicity and ΔCFS modelling, especially in regions exhibiting high static stress but limited recorded seismic activity.

6. Conclusions

This study quantifies the ΔCFS induced by the Mw 6.8 Liwa mainshock in 1994 and evaluates their correlation with the relocated seismicity along the southern segment of the Semangko Fault. By integrating high-resolution hypocenter relocation and ΔCFS modelling across multiple depths and fault geometries, several key findings have emerged:

The integration of earthquake relocation and ΔCFS modelling is effective for identifying zones where seismic hazard is modulated by static stress transfer, elevated within positive-ΔCFS lobes and potentially reduced within stress-shadow regions (ΔCFS < 0). We note that postseismic creep/viscoelastic relaxation may further decrease stresses with time; our static ΔCFS results represent the instantaneous post-mainshock state. These results highlight the critical importance of sustained seismotectonic monitoring in southern Sumatra, including the expansion of seismic station networks, integration with geodetic observations (e.g. GPS and InSAR), and the implementation of paleoseismic studies to reconstruct the region's earthquake history.

References

Alif, S. M., Fattah, E. I., & Kholil, M. (2020). Geodetic slip rate and locking depth of East Semangko Fault derived from GPS measurement. Geodesy and Geodynamics, 11(3), 222–228.https://doi.org/10.1016/j.geog.2020.04.002

Bellier, O., Sébrier, M., Pramumijoyo, S., Beaudouin, T., Harjono, H., Bahar, I., & Forni, O. (1997). Paleoseismicity and seismic hazard along the Great Sumatran Fault (Indonesia). Journal of Geodynamics, 24(1–4), 169–183.https://doi.org/10.1016/S0264-3707(96)00051-8

Cattin, R., Chamot-Rooke, N., Pubellier, M., Rabaute, A., Delescluse, M., Vigny, C., Fleitout, L., & Dubernet, P. (2009). Stress change and effective friction coefficient along the Sumatra–Andaman–Sagaing fault system after the 26 December 2004 (Mw 9.2) and the 28 March 2005 (Mw 8.7) earthquakes. Geochemistry, Geophysics, Geosystems, 10(3).https://doi.org/10.1029/2008GC002167

Daryono, M. R., Natawidjaja, D. H., & Sieh, K. (2012). Twin-surface ruptures of the March 2007 M > 6 earthquake doublet on the Sumatran Fault. Bulletin of the Seismological Society of America, 102(6), 2356–2367.https://doi.org/10.1785/0120110220

DeMets, C., Gordon, R. G., & Argus, D. F. (2010). Geologically current plate motions. Geophysical Journal International, 181(1), 1–80.https://doi.org/10.1111/j.1365-246X.2009.04491.x

Heidbach, O., Rajabi, M., Cui, X., Fuchs, K., Müller, B., Reinecker, J., Reiter, K., Tingay, M., Wenzel, F., Xie, F., Ziegler, M., Zoback, M.-L., & Zoback, M. (2018). The World Stress Map database release 2016: Crustal stress pattern across scales. Tectonophysics, 744, 484–498.https://doi.org/10.1016/j.tecto.2018.07.007

Hu, H., Zhao, D., Lin, J., & Pilia, S. (2023). A slab window beneath North Sumatra revealed by P-wave mantle tomography. Journal of Geophysical Research: Solid Earth, 128(6).https://doi.org/10.1029/2022JB025976

Hurukawa, N., Wulandari, B. R., & Kasahara, M. (2014). Earthquake history of the Sumatran Fault, Indonesia, since 1892, derived from relocation of large earthquakes. Bulletin of the Seismological Society of America, 104(4), 1750–1762.https://doi.org/10.1785/0120130201

Hutchings, S. J., & Mooney, W. D. (2021). The seismicity of Indonesia and tectonic implications. Geochemistry, Geophysics, Geosystems, 22(9).https://doi.org/10.1029/2021GC009812

Ishimaru, Y., Takada, Y., Ching, K. E., & Chang, W. L. (2025). Strong fault interaction in double-vergence structure: Lessons from the 2022 Yuli and 2022 Chihshang earthquakes, eastern Taiwan. Journal of Geophysical Research: Solid Earth, 130(7).https://doi.org/10.1029/2025JB031225

Irsyam, M., Sengara, I. W., Aldiamar, F., Widiyantoro, S., Triyoso, W., Natawidjaja, D. H., Kertapati, E., Meilano, I., Suhardjono, Suharyadi, R., & Ridwan, M. (2017). Peta Sumber dan Bahaya Gempa Indonesia Tahun 2017. PUSGEN, Kementerian PUPR, Jakarta.

Ito, T., Gunawan, E., Kimata, F., Tabei, T., Meilano, I., Agustan, …, & Sugiyanto, D. (2016). Co-seismic offsets due to two earthquakes (Mw 6.1) along the Sumatran Fault system derived from GNSS measurements. Earth, Planets and Space, 68(1), 57.

Kennett, B. L. N., Engdahl, E. R., & Buland, R. (1995). Constraints on seismic velocities in the Earth from travel times. Geophysical Journal International, 122(1), 108–124.https://doi.org/10.1111/j.1365-246X.1995.tb03540.x

King, G. C. P., & Devès, M. H. (2015). Fault interaction, earthquake stress changes, and the evolution of seismicity. In Treatise on Geophysics (2nd ed., Vol. 4, pp. 243–271). Elsevier.https://doi.org/10.1016/B978-0-444-53802-4.00077-4

King, G. C., Stein, R. S., & Lin, J. (1994). Static stress changes and triggering earthquakes. Bulletin of the Seismological Society of America, 84(3), 935–953.

Kurniawan, A., Wibowo, R. C., Nurrochman, A., Zaenudin, A., Wibowo, A., Destawan, R., & Hanif, I. (2024). Analisis seismisitas dan mekanisme fokus Sesar Sumatera: Studi kasus segmen Semangko dan Selat Sunda. Jurnal Teknologi dan Inovasi Industri (JTII), 5(2).https://doi.org/10.23960/jtii.v5i2.98

Laske, G., Masters, G., Ma, Z., & Pasyanos, M. (2013, April). Update on CRUST1.0—A 1-degree global model of Earth’s crust. Geophysical Research Abstracts, 15, 2658.

Lubis, A. M., Salman, R., Hermawan, I., Bradley, K., Feng, L., Qiu, Q., Sahputra, R., Natawidjaja, D. H., Sieh, K., & Hill, E. M. (2024). Slip rates and locking depths of the southern Sumatran Fault Zone revealed by new campaign GPS observations. Geophysical Journal International, 239(1), 248–257.https://doi.org/10.1093/gji/ggae257

Lubis, R. F., Bradley, K. E., Natawidjaja, D. H., & Bellier, O. (2024). Coseismic and postseismic deformation along the southern Great Sumatran Fault derived from InSAR analysis. Journal of Geophysical Research: Solid Earth, 129(2), e2023JB027432.https://doi.org/10.1029/2023JB027432

Mildon, Z., Toda, S., Walker, J., & Roberts, G. (2018). Evaluating models of Coulomb stress transfer: Is variable fault geometry important? Geophysical Research Letters, 43(24), 1–8.https://doi.org/10.31223/OSF.IO/7XU2B

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).https://doi.org/10.1186/s40562-022-00259-2

Natawidjaja, D. H. (2018). Updating active fault maps and slip rates along the Sumatran Fault Zone, Indonesia. IOP Conference Series: Earth and Environmental Science, 118, 012001.https://doi.org/10.1088/1755-1315/118/1/012001

Natawidjaja, D. H., & Triyoso, W. (2007). The Sumatran fault zone—from source to hazard. Journal of Earthquake and Tsunami, 1(1).https://doi.org/10.1142/S1793431107000031

Okada, Y. (1992). Internal deformation due to shear and tensile faults in a half-space. Bulletin of the Seismological Society of America, 82(2), 1018–1040.https://doi.org/10.1785/BSSA0820021018

Parija, M. P., Kumar, S., Tiwari, V. M., Biswal, S., Biswas, A., & Velliyidathu, A. (2021). Coulomb stress modeling and seismicity in the western Himalaya, India, since 1905: Implications for the incomplete ruptures of the Main Himalayan Thrust. Tectonics, 40(9), 1–16.https://doi.org/10.1029/2020TC006204

Pope, N., & Mooney, W. D. (2020). Coulomb stress models for the 2019 Ridgecrest, California earthquake sequence. Tectonophysics, 791, 228555.https://doi.org/10.1016/j.tecto.2020.228555

Pranata, B., Ramdhan, M., Hanif, M., Sulaiman, M. I., Maulana, M. P., Wandono, Widiyantoro, S., Suhardja, S. K., Hidayat, E., Supendi, P., Kusnandar, R., & Setyonegoro, W. (2023). Seismic imaging beneath Sumatra Island and its surroundings from local–regional P-wave tomography. Rudarsko-Geološko-Naftni Zbornik, 38(3), 119–132.https://doi.org/10.17794/rgn.2023.3.10

Rafie, M. T., Cummins, P. R., Sahara, D. P., Widiyantoro, S., Triyoso, W., & Nugraha, A. D. (2021). Variations in forearc stress and changes in principal stress orientations caused by the 2004–2005 megathrust earthquakes in Sumatra, Indonesia. Frontiers in Earth Science, 9, 712144.https://doi.org/10.3389/feart.2021.712144

Rafie, M. T., Sahara, D. P., Cummins, P. R., Triyoso, W., & Widiyantoro, S. (2023). Stress accumulation and earthquake activity on the Great Sumatran Fault, Indonesia. Natural Hazards, 116(3), 3401–3425.https://doi.org/10.1007/s11069-023-05816-2

Sahara, D. P., Widiyantoro, S., & Irsyam, M. (2018). Stress heterogeneity and its impact on seismicity pattern along the equatorial bifurcation zone of the Great Sumatran Fault, Indonesia. Journal of Asian Earth Sciences, 164, 1–8.https://doi.org/10.1016/j.jseaes.2018.06.002

Salman, R., Lindsey, E. O., Feng, L., Bradley, K., Wei, S., Wang, T., Daryono, M. R., & Hill, E. M. (2020). Structural controls on rupture extent of recent Sumatran Fault Zone earthquakes, Indonesia. Journal of Geophysical Research: Solid Earth, 125(2), e2019JB018101.https://doi.org/10.1029/2019JB018101

Segou, M., & Parsons, T. (2020). A new technique to calculate earthquake stress transfer and to probe the physics of aftershocks. Bulletin of the Seismological Society of America, 110(2), 863–873.https://doi.org/10.1785/0120190033

Sieh, K., & Natawidjaja, D. (2000). Neotectonics of the Sumatran Fault, Indonesia. Journal of Geophysical Research: Solid Earth, 105(B12), 28295–28326.https://doi.org/10.1029/2000JB900120

Soehaimi, A., Widarto, D. S., Matsuryono, & Effendi, I. (2002). The seismotectonic database as main parameter for prediction of the tectonic earthquake hazard level at Liwa, West Lampung District. In Proceedings of the 31st Annual Convention of the Indonesian Association of Geologists.

Sukrungsri, S., Khamsiri, S., & Pailoplee, S. (2024). Investigation of co-seismic stress and aftershock distribution along the Sumatra–Andaman subduction zone. Geoscience Letters, 11(1), 1–15.https://doi.org/10.1186/s40562-024-00366-2

Supendi, P., Rawlinson, N., Prayitno, B. S., Sianipar, D., Simanjuntak, A., Widiyantoro, S., Palgunadi, K. H., Kurniawan, A., Shiddiqi, H. A., Nugraha, A. D., Sahara, D. P., Daryono, D., Triyono, R., Adi, S. P., Karnawati, D., Daniarsyad, G., Ahadi, S., Fatchurochman, I., Anugrah, S. D., Heryandoko, N., Sudrajat, A. (2022). A previously unidentified fault revealed by the February 25, 2022 (Mw 6.1) Pasaman Earthquake, West Sumatra, Indonesia. Physics of the Earth and Planetary Interiors, 334, 106973.https://doi.org/10.1016/j.pepi.2022.106973

Supendi, P., Widiyantoro, S., Rawlinson, N., Wibowo, A., Priyobudi, P., Palgunadi, K. H., Nugraha, A. D., Imran, I., Marliyani, G. I., Daryono, D., Prayitno, B. S., Sadly, M., Karnawati, D., Sari, N., & Sugiharto, A. (2022). Analysis of the 2021 Semangko Bay earthquake sequence in southern Sumatra, Indonesia, using broadband seismic network data. Seismological Research Letters, 93(3), 1373–1381.https://doi.org/10.1785/0220210304

Tozer, B., Sandwell, D. T., Smith, W. H. F., Olson, C., Beale, J. R., & Wessel, P. (2019). Global bathymetry and topography at 15 arc sec: SRTM15+. Earth and Space Science, 6(10), 1847–1864.https://doi.org/10.1029/2019EA000658

Toda, S., Stein, R. S., Richards-Dinger, K., & Bozkurt, S. B. (2005). Forecasting the evolution of seismicity in southern California: Animations built on earthquake stress transfer. Journal of Geophysical Research: Solid Earth, 110(B5).https://doi.org/10.1029/2004JB003415

Toda, S., Stein, R. S., Sevilgen, V., & Lin, J. (2011). Coulomb 3.3—Graphic-rich deformation and stress-change software for earthquake, tectonic, and volcano research and teaching (User Guide). USGS Open-File Report 2011-1060.

Valentini, G., Volatili, T., Galli, P., & Tondi, E. (2024). Investigating the last millennium Coulomb stress transfer in the Central Apennine Fault System (CAFS). Tectonics, 43(10), e2024TC008328.https://doi.org/10.1029/2024TC008328

Waldhauser, F. (2001). hypoDD—A program to compute double-difference hypocenter locations (version 1.0-03/2001). USGS Open-File Report 01-113.https://doi.org/10.3133/ofr01113

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.https://doi.org/10.1785/0120000006

Wang, Y., Gao, Y., Morley, C. K., Seagren, E. G., Qian, X., Rimando, J. M., et al. (2023). Pleistocene accelerated exhumation within the Sumatran Fault: Implications for late Cenozoic evolution of Sumatra (Indonesia). Geophysical Research Letters, 50, e2022GL100028.https://doi.org/10.1029/2022GL100028

Wells, D. L., & Coppersmith, K. J. (1994). New empirical relationships among magnitude, rupture length, rupture width, rupture area, and surface displacement. Bulletin of the Seismological Society of America, 84(4), 974–1002.https://doi.org/10.1785/BSSA0840040974

Wessel, P., Luis, J. F., Uieda, L., Scharroo, R., Wobbe, F., Smith, W. H. F., & Tian, D. (2019). The Generic Mapping Tools version 6. Geochemistry, Geophysics, Geosystems, 20(11), 5556–5564.https://doi.org/10.1029/2019GC008515

Widiwijayanti, C., Déverchère, J., Louat, R., Sébrier, M., Harjono, H., Diament, M., & Hidayat, D. (1996). Aftershock sequence of the 1994 Mw 6.8 Liwa earthquake (Indonesia): Seismic rupture process in a volcanic arc. Geophysical Research Letters, 23(21), 3051–3054.https://doi.org/10.1029/96GL02048

Wulandari, A., Supendi, P., & Natawidjaja, D. H. (2023). Kinematic behavior of the Angkola–Sianok fault system based on the 2022 Pasaman earthquake sequence. Journal of Asian Earth Sciences, 254, 105153.https://doi.org/10.1016/j.jseaes.2023.105153

Wulandari, R., Chan, C. H., & Wibowo, A. (2023). The 2022 Mw 6.2 Pasaman, Indonesia, earthquake sequence and its implication for seismic hazard in central-west Sumatra. Geoscience Letters, 10(1).https://doi.org/10.1186/s40562-023-00279-6

Xiong, X., Shan, B., Zhou, Y. M., Wei, S. J., Li, Y. D., Wang, R. J., & Zheng, Y. (2017). Coulomb stress transfer and accumulation on the Sagaing Fault, Myanmar, over the past 110 years and its implications for seismic hazard. Geophysical Research Letters, 44(10), 4781–4789.https://doi.org/10.1002/2017GL072770

Author’s contribution

Andri Kurniawan conceived and designed the study, performed Coulomb stress modelling, interpreted results, and drafted the original manuscript. Nicolas Silaen contributed to data curated and processed seismic data (phase picks, double-difference relocation), relocation processing, and figure preparation, and assisted with writing. Khalil Ibrahim contributed to seismotectonic interpretation and manuscript editing. Arsy Nurrochman contributed to data processing, visualization, and figure preparation. Ridho Destawan contributed to seismic-catalogue preparation, quality control, and manuscript review. Nanda Hanyfa Maulida compiled datasets, performed preliminary statistical analyses, and assisted with manuscript preparation. Irfan Prasetyo contributed to parts of the modelling workflow and manuscript review. Ahmad Zaenudin verified the catalogue, contributed to data cleaning and validation, and reviewed the manuscript. Adhi Wibowo (BMKG) provided BMKG catalogue and station metadata, contributed to results validation, and reviewed the final manuscript.

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