Robust Estimation of Earthquake Magnitude in Indonesia Using PGD Scaling Law from Regional High-Rate GNSS Data

1. Introduction

Indonesia is very susceptible to earthquakes and tsunamis due to its position at the intersection of four significant tectonic plates: the Indo-Australian, Eurasian, Philippine, and Pacific [1,2,3,4]. This complex tectonic interaction creates a highly dynamic and seismically active environment, resulting in significant loss of life and economic devastation. For example, the loss of life and economic loss due to the 2024 Aceh earthquake and tsunami, it is estimated that the death toll reached 173,741 people and the financial loss was worth 4.9 billion USD [5]. The Indonesian Tsunami Early Warning System (InaTEWS) was built in 2008 to warn early of tsunamis caused by significant offshore earthquakes, lowering risks and boosting community resilience.

Earthquake magnitude estimation is crucial for earthquake and tsunami early warnings, while magnitude scaling is derived primarily from seismic sensors. However, magnitude saturation often occurs for large earthquakes (Mw > 7.5), leading to underestimated initial values. For example, in the 26 December 2004 Aceh tsunami, the Pacific Tsunami Warning Center (PTWC) initially estimated a magnitude of Mw 8.0, 11 minutes after origin time (OT), 1 hour later rise to Mw 8.5, and a few days later revised to Mw 9.2 after more seismic data became available [6,7,8]. Similarly, during the 2011 Tohoku tsunami in Japan, the Japan Meteorological Agency (JMA) first estimated Mw 7.9 within three minutes, later updating it to Mw 8.4 and finally to Mw 9.0 after several days [9,10]. If the actual size of earthquakes had been determined earlier, more effective evacuations could have been possible.

Magnitude saturation occurs when seismic signals from surrounding sensors are clipped, resulting in off-scale amplitudes. Large earthquakes generate signals that surpass the seismometer's dynamic range, resulting in incorrect results [11,12,13,14]. Furthermore, acceleration recordings are utilized in the near-field of major earthquakes because they do not undergo amplitude saturation. However, identifying the magnitude involves integrating acceleration to obtain displacement, which was caused by baseline drift due to factors such as instrument tilt and rotation [13,15,16]. High-rate GNSS data, however, provide direct displacement measurements without these limitations, making them valuable for real-time earthquake magnitude estimation [13,14,15,16].

Numerous studies have demonstrated that peak ground displacement (PGD) derived from high-frequency (1 Hz) GNSS data can accurately estimate earthquake magnitude, even before the rupture process is completed [17,18,19,20]. PGD is defined as the maximum displacement observed in a time series. Crowell et al. [12] introduced the PGD Scaling Law, which establishes a quantitative relationship between peak ground displacement and the distance from the earthquake epicenter. This algorithm enables the rapid and accurate calculation of earthquake moment magnitude (Mw) by comparing the maximum ground displacements recorded at multiple GNSS stations.

Indonesia's Tsunami Early Warning System (InaTEWS) relies on seismometer recordings to determine earthquake parameters, including location, magnitude, and moment tensor. Despite Indonesia having 473 continuously operating GNSS stations within the Indonesian Continuously Operating Reference Station (InaCORS) network, high-rate GNSS (HR-GNSS) data has not yet been incorporated into earthquake parameter determination in InaTEWS. Instead, GNSS data is primarily used in surveying and applications. [21]

This paper introduces a novel PGD Scaling Law that utilizes PGD extracted from GNSS displacement waveforms to determine earthquake magnitude in Indonesia. Initially, we examine the attenuation relationship of PGD concerning hypocentral distance using 87 displacement records from 21 moderate to large earthquakes in Indonesia. We then use this attenuation relationship to derive a PGD scaling law. Finally, we apply the PGD scaling law to estimate and compare earthquake magnitudes with the reported moment magnitudes. Our regional PGD Scaling Law, derived from Indonesia-specific GNSS data, demonstrates accuracy compared to global models, as evidenced by a more favorable mean absolute deviation (MAD) value, and provides robust estimates of earthquake magnitude.

2. Materials and Method

2.1. High-Rate GNSS Data Collection

We acquired 1 Hz high-rate GNSS data of 21 moderate to large earthquakes in Indonesia between 2009 and 2022, with magnitudes ranging from Mw 5.6 to 8.4 (Figure 1 and Table 1). This high-rate GNSS data, a significant component of our research, was collected from two networks: the Sumatran GPS Array (SuGAr) and the Indonesia Continuously Operating Reference Station (InaCORS). SuGAr stretches over a thousand kilometers at the tectonic border where the Indo-Australian and Asian plates converge. The Earth Observatory of Singapore (EOS) manages this network, which has recently extended to include 60 continuous GNSS stations. It offers essential data for understanding the Sunda megathrust and the Sumatran fault [23]. The Geospatial Information Agency (BIG) administers InaCORS, which comprises 473 GNSS stations across Indonesia. This network may enable various applications, including high-precision mapping, navigation, and geodynamic monitoring [21].

Due to the inherent noise levels of approximately 1-2 cm in GNSS data [15], stations farther from the earthquake hypocenter often fail to record substantial signals for certain moderate-magnitude events. Consequently, we selected GNSS waveforms exhibiting peak amplitudes of no less than 2 cm, resulting in 87 high-rate GNSS recordings across three components (east, north, and vertical) from 21 moderate to large earthquakes in Indonesia. The HR GNSS data, provided in the Receiver Independent Exchange Format (RINEX), were processed to derive earthquake displacement data using the Precise Point Positioning with Ambiguity Resolution (PPP-AR) technique [24,25,26], implemented through the PRIDE PPP-AR software developed by the GNSS Research Center at Wuhan University [27]. The PPP-AR method necessitates the integration of precise satellite orbit and clock corrections, thereby ensuring centimeter-level accuracy in the displacement measurements.

2.2. PGD Scaling Law

Peak Ground Displacement (PGD) is the maximum dynamic displacement recorded at a GNSS station during an earthquake. This value is calculated as the Euclidean norm of the three components (N, E, U), expressed mathematically as follows:

$$\mathrm{P}\mathrm{G}\mathrm{D}=\mathrm{m}\mathrm{a}\mathrm{x}\left(\sqrt{{\mathrm{N}\left(\mathrm{t}\right)}^2+{\mathrm{E}\left(\mathrm{t}\right)}^2+{\mathrm{U}\left(\mathrm{t}\right)}^2}\right)$$

(1)

PGD Scaling Law defines the relation between PGD, moment magnitude (Mw), and the distance from the earthquake hypocenter to the GNSS station (R) as follows [15,16,17,18].

$$\mathrm{l}\mathrm{o}\mathrm{g}\mathrm{P}\mathrm{G}\mathrm{D}=\mathrm{A}+\mathrm{B}.\mathrm{M}\mathrm{w}+\mathrm{C}.\mathrm{M}\mathrm{w}.\mathrm{l}\mathrm{o}\mathrm{g}\left(\mathrm{R}\right)$$

(2)

$$\mathrm{M}\mathrm{w}={\displaystyle \frac{\mathrm{l}\mathrm{o}\mathrm{g}\mathrm{P}\mathrm{G}\mathrm{D}-\mathrm{A}}{\mathrm{B}+\mathrm{C}.\mathrm{l}\mathrm{o}\mathrm{g}\left(\mathrm{R}\right)}}$$

(3)

where M_w is moment magnitude (M_w) determined by BMKG; R is the hypocentral distance in km; and A, B and C are the regression coefficients.

Another study previously employed a PGD Scaling Law for moderate to large earthquakes, based on global displacement HR GNSS worldwide [13,14,15]. We attempted to develop our local PGD scaling law using 87 HR Data GNSS from 21 moderate to major earthquakes in Indonesia. We use linear regression with the least squares method to analyze the 87 PGD measurement data and obtain the regression coefficients (A, B, C). We estimate them using a bootstrap approach, randomly removing 10% of PGD measurements and rerunning the regression. We repeat this step 1000 times to estimate the variance of the coefficients and present the uncertainties as the 95% confidence intervals of each coefficient.

3. Results

To evaluate the effectiveness of early peak ground displacement (PGD) measurements for rapid magnitude estimation, the kinematic displacement time series was truncated to a 5–7 minute window following the earthquake origin time (OT). Figure 2 presents the displacement waveforms of the Mw 7.7 Sinabang earthquake (6 April 2010) in three components (East, North, and Up), as recorded by six high-rate GNSS stations. A systematic decrease in displacement amplitude with increasing hypocentral distance is observed, consistent with theoretical expectations of ground motion attenuation. Figures S1–S5 in the Supplementary Material shows the epicenters and GNSS station locations alongside the displacement waveforms recorded by at least six High-Rate GNSS stations for selected earthquakes. Table S1 in the Supplementary Material presents the Peak Ground Displacement (PGD) calculation results and the corresponding hypocentral distances for each HR GNSS station across 21 moderate to large earthquakes. The PGD values, derived from 87 HR GNSS stations, range from 2.01 cm to 292.96 cm. The smallest PGD (2.01 cm) was recorded at the CPMK station during the Mw 6.9 earthquake on December 15, 2017, in Tasikmalaya, West Java. In contrast, the largest PGD (292.96 cm) was observed at the PALP station for the Mw 7.5 earthquake on September 28, 2018, in Palu-Donggala, Central Sulawesi. The hypocentral distances of the HR GNSS stations ranged from 17 km to 1287 km.

The regression coefficients obtained from the linear regression analysis were as follows: A = -4.729 ± 0.255, B = 1.055 ± 0.045, and C = -0.121 ± 0.006, with a residual standard deviation of 0.30 magnitude units. The Peak Ground Displacement (PGD) Scaling Law was derived based on data from 87 High-Rate GNSS stations, corresponding to 21 moderate to major earthquakes in Indonesia, and is expressed as follows:

$$\mathrm{l}\mathrm{o}\mathrm{g}\mathrm{P}\mathrm{G}\mathrm{D}=-4.729+1.055.\mathrm{M}\mathrm{w}+-0.121.\mathrm{M}\mathrm{w}.\mathrm{l}\mathrm{o}\mathrm{g}\left(\mathrm{R}\right)$$

(4)

Figure 3 presents the distribution of PGD as a function of hypocentral distance. PGD measurements at each station are plotted, with color and marker types distinguishing different earthquakes. The oblique lines illustrate predicted magnitudes as a function of PGD and hypocentral distance. The PGD measurements for each event generally fluctuate among the oblique lines and decrease in amplitude as the hypocentral distance increases. The plot reveals a clear attenuation trend, with PGD decreasing approximately logarithmically with distance. Events of higher magnitude consistently produce larger PGDs at all distances.

Several well-recorded large-magnitude earthquakes (e.g., the doublet Simeulue 2012, Padang 2009, Sinabang 2010, Palu-Donggala 2018) form an upper bound to the PGD distribution. Meanwhile, moderate-magnitude events (e.g., Cianjur 2022, Mw 5.6) show PGDs confined to the lower part of the graph. Additionally, shallow earthquakes, particularly strike-slip events such as Palu 2018 (Mw 7.5), generate high PGDs at relatively short hypocentral distances, highlighting the hazard potential of near-field ruptures. In contrast, deeper events (e.g., Tasikmalaya 2017, depth 107 km) generate significantly lower surface displacements, even at moderate distances.

Using the determined PGD Scaling Law, we compute the mean PGD magnitude (M_PGD) estimation and its standard deviation for each of the 21 events, as detailed in Table 1. Table 1 also informs 21 earthquake events across Indonesia between 2009 and 2022 that were analyzed to evaluate the relationship between moment magnitude (Mw) and magnitude estimated from peak ground displacement (M_PGD). The events exhibit diverse tectonic settings, encompassing strike-slip, reverse, and normal fault mechanisms. Most earthquakes occurred at shallow depths (<30 km), with a few exceptions such as the Padang (81 km) and Tasikmalaya (107 km) events, which occurred at intermediate depths.

Table 1 also shows the number of High-Rate GNSS (HR GNSS) stations that recorded surface displacement during earthquakes, ranging from 1 to 19 per event. As of the most recent data, Indonesia has 473 continuous HR GNSS stations under the InaCORS (Indonesian Continuous GNSS Network), although the overall density remains relatively sparse. The average inter-station distance is approximately 72 km, and has steadily expanded since 2021 [29], while adequate for broad-scale monitoring, may limit the spatial resolution of surface displacement data, particularly in seismically active regions.

The number of GNSS stations significantly influenced the robustness of the M_PGD estimates. Events with more than five GNSS stations, such as the Simeulue and Mentawai earthquakes, exhibited lower standard deviation in M_PGD. In contrast, events monitored by only one or two stations, like the Sangihe or Palu 2017 earthquakes, showed larger standard deviation ranges in magnitude estimation. Overall, M_PGD estimates are closely aligned with Mw, with a mean residual standard deviation of 0.3 magnitude units, indicating the robustness of PGD scaling for rapid magnitude estimation.

We also conducted retrospective magnitude evolution estimations for several large earthquakes using high-rate GNSS data (Figure 4), focusing on events recorded by a minimum of six stations to ensure magnitude stability, with a mean error of ≤ 0.3 magnitude units as found in research by Ruhl et al. [17]. The five earthquakes include: the Mw 7.7 Sinabang (2010), Mw 7.2 Mentawai (2010), Mw 8.4 Simeulue (2012), Mw 8.1 Simeulue (2012), and Mw 7.5 Palu-Donggala (2018) events.

M_PGD estimates initially exhibit significant variability within the first 50–100 seconds from OT for all events due to the dynamic rupture evolution and limited early PGD observations. As more GNSS data become available and peak displacements are recorded, the magnitude estimates converge progressively toward stable values. In the case of the Sinabang (2010) and Mentawai (2010) earthquakes, the M_PGD stabilized within approximately 125 s from OT, accurately matching the final Mw values of 7.7 and 7.1, respectively. Similarly, for the doublet Simeulue earthquakes (2012), the stable M_PGD was achieved at 8.53 ± 0.15 and 8.03 ± 0.12 within approximately 300 s from OT, which closely correspond to the final moment magnitude from BMKG Mw 8.4 and Mw 8.1, respectively. Across all events, BMKG initial seismic magnitude estimates (red stars) tend to underestimate the final Mw, especially for the great magnitude events (Mw > 8.0). In contrast, PGD-based estimates show a more consistent trend toward the final values. These results demonstrate the robustness and reliability of the PGD scaling method for rapid earthquake magnitude estimation.

4. Discussion

Based on the PGD (Peak Ground Displacement) magnitude estimation results for the 21 earthquakes used to develop the PGD scaling law, the standard deviation of the estimated magnitudes ranges from 0.02 to 0.74 (Table 1). Generally, lower standard deviation values are observed for earthquakes in which at least six GNSS stations recorded the displacement. For instance, significant events along the Sumatra plate boundary, such as the Mw 7.7 Sinabang earthquake (2010), the Mw 8.4 Simeulue earthquake (2012), and the Mw 8.1 Simeulue earthquake (2012), exhibited more minor uncertainties. This finding is consistent with the results reported by Ruhl et al. [17], who demonstrated that stable magnitude estimates, with uncertainties equal to or less than 0.3 magnitude units, require observations from a minimum of six high-rate GNSS stations. It also aligns with the finding by Gao et al. [16], who noted that an increase in GNSS stations improves the accuracy of the final magnitude estimates. Notably, the estimations for the two great earthquakes in Simeulue (Mw 8.4 and Mw 8.1), with M_PGD values of 8.53 ± 0.15 and 8.03 ± 0.12, respectively, confirm that the PGD magnitude estimations do not exhibit saturation even for very large events.

Anomalously high PGD magnitude estimates were observed for the Mw 7.5 Palu earthquake (28 September 2018), reaching 7.98 ± 0.74. This overestimation is primarily attributed to the large peak ground displacement (2.93 m; Table S1) recorded at the PALP GNSS station, located approximately 78 km south of the epicenter. The rupture propagated southward along the Palu-Koro fault at supershear velocity, exceeding the crust's shear wave speed, resulting in intense fault-parallel ground motion [30,31]. PALP, situated near the fault trace, recorded a fault-parallel velocity (~1.0 m/s) substantially greater than the fault-normal component (~0.7 m/s) [32], consistent with theoretical and numerical models of supershear rupture dynamics [33]. The location of the epicenter and PALP station, along with the observed displacement, are shown in Figure 5. The retrospective analysis reveals that the mean M_PGD reached 7.5 at 78 s after OT and continued to increase, eventually stabilizing at M_PGD 7.98 around 126 s from OT. This finding highlights the capability of the PGD method to provide a rapid and stable magnitude rather than seismic magnitude estimates, which remained underestimated during the initial time of the earthquake.

Numerous studies have established Peak Ground Displacement (PGD) Scaling Laws utilizing high-rate GNSS (HR-GNSS) data from global seismic events, including those by Melgar et al. [13], Crowell et al. [14], and Ruhl et al. [15]. Building upon this foundation, the present study aims to develop a PGD Scaling Law based on regional HR-GNSS displacement data specifically recorded from seismic events in Indonesia. Table 2 compares the PGD Scaling Law coefficients derived in this study and those reported in the previous works.

To evaluate the accuracy of the proposed PGD Scaling Law in estimating the magnitude from peak ground displacement (MPGD), a comparative analysis was conducted using 21 moderate to large earthquakes that occurred in Indonesia. The mean absolute deviation (MAD) between the MPGD-derived magnitudes and the catalog moment magnitudes (Mw) was calculated for this study and for scaling laws proposed by Melgar et al. [13], Crowell et al. [14], and Ruhl et al. [15]. As shown in Figure 6, the MAD for the proposed model is 0.21, which is lower than the MADs of 0.29, 0.39, and 0.34 for the models by Melgar et al., Crowell et al., and Ruhl et al., respectively.

Based on Indonesia-specific high-rate GNSS data, the regional PGD Scaling Law developed in this study demonstrates improved accuracy over global models, as reflected by a lower mean absolute deviation (MAD), and yields robust and reliable earthquake magnitude estimates. These results indicate that the PGD Scaling Law derived from this study produces a more accurate fit that captures the regional tectonic characteristics of the Indonesian region.

Based on the standard operating procedure of the Indonesia Tsunami Early Warning System (InaTEWS), since 2024, real-time earthquake parameters such as origin time, hypocenter location, and magnitude must be provided by BMKG in less than 3 minutes [34]. From the results of retrospective analysis of several large earthquakes, the magnitude estimation results were obtained with stable GNSS high rate PGD data in less than 3 minutes, for example the Sinabang earthquake, Mentawai earthquake, and Palu-Donggala earthquake, while for the Simeleu doublet earthquake where the hypocenter location in ocean quite far from the Sumatra mainland where the GNSS sensor is located, around 300 km so that a stable magnitude estimate can be obtained approximately 300 seconds after OT. The results reaffirm the capability of GNSS displacement data to complement and enhance traditional seismic networks to provide reliable earthquake magnitude estimates that are crucial in tsunami and earthquake early warning systems, particularly for large and tsunamigenic events in the Indonesian region. This proves that PGD Magnitude has great potential to be implemented in real time in the Indonesian Tsunami Early Warning System. However, challenges such as station density and data latency remain key considerations and future works for real-time operational implementation.

5. Conclusions

We propose a novel regionally calibrated Peak Ground Displacement (PGD) Scaling Law for rapid and reliable earthquake magnitude estimation in Indonesia. This model is developed using 87 high-rate GNSS (HR-GNSS) displacement records from 21 moderate to large earthquakes and establishes a robust empirical relationship between PGD, hypocentral distance, and moment magnitude (Mw). The resulting PGD-derived magnitudes (MPGD) demonstrate strong agreement with reported catalog Mw values, achieving a mean absolute deviation (MAD) of 0.21, outperforming previously published global models.

Notably, the proposed scaling law shows no evidence of magnitude saturation, even for great earthquakes (Mw > 8.0), and maintains accuracy across various tectonic regimes, including strike-slip, reverse, and normal fault mechanisms. Retrospective analyses further confirm that for events recorded by six or more GNSS stations, PGD-based magnitude estimates converge rapidly, typically within 2–3 minutes of origin time, meeting the operational requirements of Indonesia’s Tsunami Early Warning System (InaTEWS). These findings reaffirm the value of HR-GNSS displacement data as a complementary tool to seismic networks for real-time magnitude estimation, particularly in tectonically complex and tsunami-prone regions like Indonesia. For continuous improvement, further validation using real-time HR GNSS displacement data from future seismic events in Indonesia will be essential to evaluate the scalability and operational reliability of the proposed model.