The previous studies of the Sumatra region have shown that the Bengkulu region is a highly active seismic zone due to the collision of the Indo-Australian and Eurasian tectonic plates (e.g., [Chlieh et al., 2008; Prawirodirdjo et al., 2010]). The Australian plate is converging with the southeastern segment of the Eurasian Plate (Sundaland block) at a plate motion rate of ~59 mm/yr [Bock et al., 2003]. The 2020 Bengkulu earthquake occurred within the subduction zone between these two plates, where the Indo-Australian Plate is subducting beneath the Eurasian Plate.

The Bengkulu region, located in the western coast of Sumatra, Indonesia, has experienced several significant earthquakes in the past (e.g., [Natawidjaja et al., 2006]). Historically, many large earthquakes have occurred in central Sumatra; the Mentawai segment produced MW 8.7–8.9 and 8.9–9.1 events in 1797 and 1833, respectively [Newcomb, McCann, 1987; Natawidjaja et al., 2006]. One of the most notable events was the 2007 Bengkulu earthquake, which had a magnitude of MW 8.4 and caused significant damage to infrastructure and loss of life [Gusman et al., 2010; Konca et al., 2008; Lubis et al., 2013; Zheng et al., 2018]. Then, the MW 7.8 Mentawai earthquake of October 25, 2010 caused a tsunami [Hill et al., 2012; Li, Huang, 2013; Tsang et al., 2016]. In addition to coseismic deformation, this region is also affected by ongoing postseismic deformation resulting from the 2004 Sumatra–Andaman earthquake [Alif et al., 2024]. Another earthquake with a magnitude of MW 6.8 occurred on August 18, 2020 at 22:23:59 UTC in the Bengkulu region. The United States Geological Survey (USGS) reported that the 2020 Bengkulu earthquake epicenter was located approximately ~130 km southwest of the city of Bengkulu (4.322° S; 101.135° E) at a depth of 22 km depth. The focal mechanisms of the earthquake suggest two possible directions of fault source (strike 313° and 127°) (https://earthquake.usgs.gov/earthquakes/search/). However, the faults of the 2020 Bengkulu earthquake are poorly understood.

GNSS is a very reliable tool for measuring the dynamics of the earth’s crust [Heliani et al., 2024; Ansari et al., 2024], including volcanic deformation [Anggara et al., 2023], earthquake deformation estimation [Tanaka et al., 2019; Xiang et al., 2024], and fault slip rate calculations [Tin et al., 2022]. GNSS-based geodynamic approaches have been widely applied for different geodynamic settings, such as shear faults (San Andreas faults) [Xu et al., 2021], subduction zones in Japan [Okada, Nishimura, 2025], rift regions (Red Sea, East Africa) [Stamps et al., 2014] and coseismic effects of the Hovsgol earthquake, Mongolia [Lukhnev et al., 2022]. In this study, we aim to provide a comprehensive characterization of the coseismic slip distribution of the MW 6.8 Bengkulu earthquake by inverting continuous GNSS observations from the Indonesian Continuously Operating Reference Stations (Ina-CORS) and Sumatran GPS Array (SuGAR) networks. We specifically focus on testing two plausible fault geometries to better constrain the fault rupture parameters. Furthermore, we assess the Coulomb stress transfer resulting from the mainshock and compare it with the spatial distribution of aftershocks reported by the USGS to explore the potential triggering mechanisms. Through these objectives, the study contributes to a deeper understanding of seismic source processes and stress interactions in the Bengkulu region.

This study used the GNSS data from various sites in Bengkulu, Indonesia that were obtained from 14 stations of the Geospatial Information Agency of Indonesia (BIG) and 10 stations of the Sumatran GPS Array (SuGAR) (Table 1). The GNSS sites were located around the Mentawai island and the coast of Bengkulu (Fig. 1). The GNSS data have proven useful in estimating coseismic slip source modelling (e.g., [Alif et al., 2021; Gunawan et al., 2022]). The GNSS data were recorded every 30 seconds from 10 days before the August 18, 2020 MW 6.8 Bengkulu earthquake and 10 days after the earthquake.

The daily solution was estimated using BERNESE version 5.2 from the Astronomical Institute of the University of Bern (AIUB) [Dach et al., 2015]. To minimize errors in GNSS, use has been made of the precise satellite orbits (PRE) data, ocean loading and atmospheric loading corrections, sub-daily pole, nutation model, precise orbit file (SP3), earth rotation pole (IEP), and satellite and receiver clock corrections. The daily GNSS data were processed using the International GNSS Service (IGS) data from 10 stations (DARW, DGAR, KARR, HYDE, IISC, NTUS, PERT, PIMO, YARR, YAR2) (e.g., [Alif et al., 2025]) in the International Terrestrial Reference Frame (ITRF) 2014 [Altamimi et al., 2016]. Three components were used to calculate the coseismic deformation obtained by averaging the data for 10 days before and 10 days after the earthquake, and a coseismic slip model for the MW 6.8 Bengkulu earthquake.

The coseismic slip model was constructed using seismic offsets observed approximately 130 km from the epicenter of the earthquake. Model 1 and Model 2 were calculated based on the mechanism of the earthquake, possibly caused by faulting, using Nodal Plane 1 (NP1) and Nodal Plane 2 (NP2). Model 1 uses a strike angle (θ) of 313°, a rake angle (δ) of 70–130°, and a dip angle (λ) of 8°, a length (L) of 60 km, a width (W) of 50 km, and a patch size of 5 km. Meanwhile, Model 2 uses a strike angle (θ) of 127°, a dip angle (λ) of 82°, a length (L) of 65 km, a width (W) of 40 km, and a patch size of 5 km. We investigated two fault scenarios from the USGS thrust fault plane mechanism. We calculated the coseismic slip distribution of the earthquake using the steepest descent inversion method [Wang et al., 2009; Wang et al., 2013]:

F(s)=||Gs–y||2+α2||Hτ||2

where G – Green’s function calculated using an elastic half-space dislocation model [Okada, 1992]; s – coseismic slip on each patch, y is the apparent displacement; α2 – smoothing factor; H – Laplacian operator, and τ is the shear stress drop.

The coseismic displacement of the August 18, 2020 MW 6.8 earthquake showed that points near the station experienced greater displacement than those farther from the epicenter. The horizontal displacement values at the station near the epicenter, MLKN, were –3.66 mm for dE and –0.96 mm for dN, while at station BSAT they were 2.95 mm for dE and –1.27 mm for dN (Fig. 2, a). However, at father stations, such as CLHT, CPBM, CMTB, and CSDH, the displacement values were relatively small, of ~0.5 mm. Since CSDH, CMTB, CPBM, and CLHT stations are far from the epicenter, our coseismic model involves only closer points, directly affected by the earthquake displacement. Judging from the displacement values, the influence of displacement decreased with distance from the epicenter. The vertical displacement at station MLKN was 5.40 mm, and at station BSAT – –6.34 mm; the displacement values at vertical stations varied greatly (Fig. 2, b).

Model 1 and Model 2 are relative to the direction of subduction the thrust fault (Fig. 3). The inversion results of Model 1 yielded a seismic moment of 1.73·1019 N·m equivalent to MW 6.76 with a misfit of ~0.4 mm. However the inversion results of Model 2 yielded a seismic moment of 1.13·1019 N·m equivalent to MW 6.64 with a misfit of ~0.6 mm. The misfit between the coseismic displacement data and the model is calculated using the mean absolute error (MAE) [Gunawan et al., 2022]. The slip distribution of Model 1 is concentrated in the southern part of the earthquake epicenter; the Model 2 slip distribution pattern, also located there, has a significant dip of 82°. Both models indicate that there is the Indo-Australian plate’s interaction with the Eurasian plate in the subdiuction zone. These results are based on an assumed rigidity of 40 GPA [Gusman et al., 2010]. According to the USGS release, the seismic moment value is 1.80·1019 N·m or equivalent to MW 6.77. Among the two earthquake models tested, Model 1 demonstrates the best agreement with the USGS solution, yielding a mean absolute error (MAE) of approximately 0.4 mm in the strike direction. This model aligns with the orientation of the subduction-related NP1.

The Coulomb failure stress (∆CFS) method can be used to investigate the relationships between earthquakes that may trigger subsequent earthquakes such as mainshock-aftershock sequences, and those between tectonic and volcanic earthquakes [Zhang et al., 2012]. In this study, the MW 6.8 Bengkulu earthquake was followed by a MW 6.9 aftershock, according to USGS reports, which was then followed by aftershocks with a magnitude less than 5. Static stress triggering was calculated based on the change in Coulomb failure stress (∆CFS) [Toda et al., 2011], using a friction coefficient of 0.4 [Ishibe et al., 2017; King et al., 1994]. Based on the coseismic slip model, Model 1 shows the best agreement with the USGS release with a minimum misfit, and this was used to investigate the stress distribution. The Coulomb stress results showed that there is a relationship between the mainshock and the aftershocks, as the distribution of the aftershocks occurred in areas of increased stress, including the MW 6.9 earthquake (Fig. 4).

Based on the GPS measurements, the coseismic displacement of the August 18, 2020 earthquake was insignificant, less than ~5 mm. The best-fit model is Model 1 with a misfit value of ~0.4 mm, resulting in a seismic moment of 1.73·1019 Nm equivalent to MW 6.76, which is in agreement with the reported USGS results, suggesting a thrusting fault due to the subduction of the Indo-Australian plate. The analysis of broad-band regional seismic-displacement waveforms could highly contribute to the further development of the seismic slip source model in the Bengkulu region, Indonesia.