1 Introduction

The subduction of the Indo-Australian plate beneath the Eurasian plate off the west coast of Sumatra has often resulted in large tsunamigenic earthquakes. Great earthquakes have occurred along this plate boundary, including the 2004 Sumatra–Andaman earthquake (Mw 9.1) and the 2005 Nias earthquake (Mw 8.6).

On September 12, 2007, a great earthquake (Mw 8.4) occurred offshore the Bengkulu province of Sumatra. According to the US Geological Survey, the mainshock at 18:10 local time (11:10 UTC) was followed by two large aftershocks (Mw 7.9 and Mw 7.1) at 23:49 UTC and 03:35 UTC, respectively. The mainshock generated a large tsunami along the Bengkulu coastline and caused damage at several locations in the Bengkulu coastal area (Borrero et al. 2009). A few years later, another large earthquake (Mw 7.8) occurred on October 25, 2010, off the west coast of Mentawai Islands at 21:42 local time (14:42 UTC). This earthquake also generated a large tsunami and caused more than 400 casualties around the Mentawai Islands. Seismological analysis indicates that the 2010 Mentawai event was a “tsunami earthquake” which had an unusually long duration of 100–120 s (Lay et al. 2011) and produced a much larger tsunami than expected given its seismic magnitude (Satake et al. 2013). The tectonic settings and source areas of the earthquakes are shown in Fig. 1.

A previous study (Gusman et al. 2010) showed that the tsunami heights surveyed along the Central Sumatra coast resulting from the 2007 Bengkulu earthquake (Borrero et al. 2009) were well explained by the tsunami simulation based on the source model estimated using joint inversion of tsunami waveforms and geodetic data. On the other hand, the tsunami heights along the coast of Mentawai Islands resulting from the 2010 Mentawai event were clearly underestimated by the tsunami simulation based on the source model estimated using tsunami waveform inversion (Satake et al. 2013). This indicates that the tsunami heights along a coastline resulting from tsunami earthquakes such as the 2010 Mentawai tsunami earthquake are difficult to forecast. This has been an important issue in Indonesia since the 2010 Mentawai tsunami earthquake.

Tanioka et al. (2017) recently developed a method combining the W-phase inversion and scaling relationship with a depth-dependent rigidity to estimate an appropriate source model for large earthquakes, including tsunami earthquakes, in Central America. They showed that the tsunami heights observed along the Nicaraguan coast resulting from four large earthquakes, including the 1992 Nicaragua tsunami earthquake, were well explained by the computed tsunamis from the rapidly estimated fault models. The 2010 Mentawai tsunami earthquake that occurred off the west coast of Sumatra had similar characteristics as those of the 1992 Nicaragua tsunami earthquake.

In this study, we estimate appropriate rigidities, such as the depth-dependent rigidity used by Tanioka et al. (2017), to determine appropriate source models for tsunami height estimations of the 2007 Bengkulu earthquake (Mw 8.4) and the 2010 Mentawai tsunami earthquake (Mw 7.8). Additionally, we determine the depth-dependent rigidity relation for Central Sumatra in our tsunami height forecasting method.

2 Data

The broadband seismic waveforms used in this study were downloaded from the Incorporated Research Institutions for Seismology (IRIS) Data Management Center (https://ds.iris.edu/ds/nodes/dmc/data/). The vertical components of the seismic waveforms observed at the stations at an epicentral distance between 5° and 50° were used to perform W-phase inversion.

The bathymetric dataset from the General Bathymetric Chart of the Oceans (GEBCO; https://www.gebco.net/data_and_products/gridded_bathymetry_data/), with a 30″ resolution was used for the tsunami simulation. The topographic data was downloaded from the Seamless Digital Elevation Model (DEMNAS; tides.big.go.id/DEMNAS/) for Mentawai Islands. The DEMNAS data were then manually corrected to create realistic topographic data with an interpolation if the observed tsunami heights of the 2010 Mentawai tsunami (Satake et al. 2013) were less than the DEMNAS elevation data at the observed points. The corrected elevation from the coast to the observed points then became less than the observed tsunami heights. The coastline was also checked against the Google Maps platform.

Less than 3 days after the 2007 Bengkulu earthquake, Borrero et al. (2009) conducted a post-tsunami field survey that covered 300 km of the Sumatra coast centered on the city of Bengkulu. For the 2010 Mentawai tsunami earthquake, Satake et al. (2013) conducted a post-tsunami survey along the western coast of northern and southern Pagai Island of the Mentawai Islands. For both events, measured tsunami data were used for comparison with computed values. For both survey datasets, in the case of multiple measurements at the same location, the largest measurement was used in the comparison.

3 Method

3.1 W-phase Inversion

The W-phase is a distinct long-period (up to 1000-s) phase, arriving before the S phase. Because of its very long period, the W-phase is better able to represent the tsunami-generating potential of an earthquake. Benavente et al. (2016) recently developed a rapid slip inversion method using the W-phase. In this paper, the W-phase inversion algorithm (Kanamori and Rivera 2008) was used to estimate the centroid location, magnitude, and earthquake mechanism for two large earthquakes: the 2007 Bengkulu earthquake (Mw 8.4) and the 2010 Mentawai earthquake (Mw 7.8). Time-domain deconvolution was used to retrieve ground displacement data from broadband seismogram network stations, which was then bandpass-filtered between 0.001 and 0.005 Hz using a fourth-order Butterworth filter. The time window (in seconds) was set by the P arrival time + 15 times the distance (degrees) between the epicenter and station. A detailed description of the W-phase inversion used in this study is provided in Kanamori and Rivera (2008) and Duputel et al. (2011, 2012).

3.2 Fault Models for Tsunami Simulations

Tanioka et al. (2017) used a simple rectangular fault model for tsunami computation to forecast tsunami heights along the coast, as it remains difficult to estimate the slip distribution of an earthquake occurring in Central America in real time. Because similarly sparse station coverage characterizes the coast of Sumatra, we seek to establish whether a simple rectangular fault model can be effectively used for Sumatra megathrust earthquakes as well. The fault length, L (km), and fault width, W (km), were calculated from the magnitude (Mw) estimated by the W-phase inversion using the scaling relationships determined by Blaser et al. (2010), i.e., \(\log L = - \;2.28 + 0.55M_{{\text{w}}}\) and \(\log W = - \;1.8 + 0.45M_{{\text{w}}}\), respectively. The fault center was at the centroid of the earthquake estimated by the W-phase inversion. The estimated earthquake parameters of strike, dip, and rake angles from the W-phase inversion were also used.

3.3 Numerical Tsunami Simulations

The vertical ocean bottom deformation resulting from the estimated fault model was calculated using the equations of Okada (1985). This vertical ocean bottom deformation was used as the initial sea surface deformation. The tsunami propagation was then numerically computed using the nonlinear shallow-water equations as governing equations, which were solved by a finite different scheme (Goto et al. 1997). Tsunami inundation was numerically computed using a moving boundary condition (Imamura 1996). A Manning’s roughness coefficient of 0.025 m−1/3 was assumed and used in the nonlinear shallow-water equations (Gusman et al. 2014). Four nested grid systems with grid sizes of 54″, 18″, 6″, and 2″ were used to simulate the tsunami.

Fig. 1
figure 1

Tectonic settings and source areas of recent great underthrust earthquakes off the island of Sumatra: 2004 Sumatra–Andaman earthquake, 2005 Nias earthquake, 2007 Bengkulu earthquake, and 2010 Mentawai tsunami earthquake. Stars show the earthquake epicenters

Full size image

3.4 Slip Amount (Rigidity) Estimation

In general, the slip amount, D (m), of a rectangular fault model is calculated from the following equations (Kanamori and Anderson 1975; Hanks and Kanamori 1979):

$$D = \frac{{M_{0} }}{\mu LW},\quad M_{0} = 10^{{1.5M_{{\text{w}}} + 9.1}}$$
(1)

where \(M_{0}\) (Nm) is the seismic moment from the result of the W-phase inversion, L is the fault length (m), W is the fault width (m), and \(\mu\) is the rigidity (N/\({\text{m}}^{2}\)). Tanioka et al. (2017) assumed depth-dependent rigidity based on the results of Bilek and Lay (1999) and calculated the slip amount using Eq. (1). In this study, we first estimated a slip amount by comparing the surveyed and computed tsunami heights. Finally, the corresponding rigidity is calculated from the equation \(\mu = \frac{{M_{0} }}{DLW}\), to use for the tsunami forecast.

To determine the best slip amount, we used the Aida number, K (Aida 1978), the geometric mean ratio between the observed and computed tsunami heights, as expressed in the following equation:

$$Log K = \frac{1}{N}\sum\nolimits_{i = 1}^{N} {\log K_{i} } ;\quad K_{i} = \frac{{O_{i} }}{{S_{i} }}$$
(2)

where \(K_{i}\) is the ratio between the observed tsunami heights (\(O_{i}\)) and the simulated tsunami heights (\(S_{i} )\) at each location i. If K is less than 1.0, the computed tsunami heights overestimate those observed. If K is greater than 1.0, the computed tsunami heights underestimate those observed. The Aida number, K, is commonly used to evaluate forecast accuracy (Gusman et al. 2014). In this study, the best slip amount is determined as that with a K value closest to 1.

4 Results

4.1 Fault Model for the 2007 Bengkulu Earthquake (Mw 8.4)

The vertical components of the broadband seismic data from 27 stations are used for W-phase inversion (Fig. 2a). The solution of the W-phase inversion provides a seismic moment of 4.34 × 1021 Nm (Mw 8.4), a thrust mechanism (strike = 332°, dip = 14°, and rake = 107°), a centroid location (latitude 4.38° S and longitude 100.99° E), and a centroid depth of 30 km. The vertical components of the observed and computed waveforms are compared in Fig. 2a. The fault length and width estimated from the scaling relationships (Blaser et al. 2010) are 208 km and 91 km, respectively (Fig. 2b). This seismic moment for the 2007 Bengkulu earthquake is consistent with that of 4.2 × 1021 Nm (Mw 8.4) estimated from tsunami waveform data analysis (Lorito et al. 2008). The fault model location is also consistent with the slip distributions estimated by Lorito et al. (2008) and Gusman et al. (2010).

Fig. 2
figure 2

Fault model estimated for the 2007 Bengkulu earthquake: a station distributions (red triangles) used in the W-phase inversion, b concatenated observed (black trace) and synthetic (red trace) W-phase waveforms for vertical component, and c mechanism estimated by the W-phase inversion and the fault model estimated from the scaling relationships. The yellow star represents the centroid location of the 2007 Bengkulu earthquake

Full size image

4.2 Tsunami Simulation for the 2007 Bengkulu Earthquake

The numerical simulations for the 2007 Bengkulu tsunami were conducted using the estimated fault model with a length of 208 km and a width of 91 km as shown in Fig. 2. For the tsunami simulation, the finest grid spacing of 2″ was used at 12 locations along the Bengkulu coastal area, as shown in Fig. 3, where tsunami heights between 1.2 and 3.5 m were surveyed by Borrero et al. (2009). The tsunami inundation simulation was not used in this computation because most of the surveyed tsunami heights were less than 3 m. The previous study by Gusman et al. (2010) also shows that the surveyed tsunami heights were well explained by the tsunami computation without the inundation computation. The tsunami heights at the coast were compared with the surveyed data. The tsunamis were computed using four different slip amounts (or corresponding rigidities for the seismic moment of 4.34 × 1021 Nm): 8.76 m (μ = 2.6 × 1010 Nm−2), 8.13 m (μ = 2.8 × 1010 Nm−2), 7.59 m (μ = 3.0 × 1010 Nm−2), and 7.12 m (μ = 3.2 × 1010 Nm−2).

Fig. 3
figure 3

Comparison of measured tsunami heights (black bar) by Borrero et al. (2009) and computed tsunami heights (red bar) using the fault model with the best slip amount of 7.59 m (corresponding rigidity of 3.0 × 1010 Nm−2) for the 2007 Bengkulu earthquake

Full size image

Figure 4 shows the Aida number, K, the geometric mean ratio between the observed tsunami heights and those computed (Eq. 2), with the slip amounts in the x-axis. The best slip amount of 7.59 m (μ = 3.0 × 1010 Nm−2) was found for the K number of 1.05. The computed tsunami heights from the best slip amount were compared to the measured tsunami heights by Borrero et al. (2009) as shown in Fig. 3 and Table 1. The measured tsunami heights between 1.2 and 3.5 m by Borrero et al. (2009) were well explained by computed tsunami heights between 0.6 and 4.0 m.

Fig. 4
figure 4

Aida number K, the geometric mean ratio between the observed and computed tsunami heights, with various slip amounts of the estimated fault model for the 2007 Bengkulu earthquake

Full size image
Table 1 Comparison of the measured and computed tsunami heights for the 2007 Bengkulu earthquake
Full size table

In this study, we found that the best source model with a slip amount of 7.59 m (μ = 3.0 × 1010 Nm−2) for the Bengkulu earthquake reproduced the observed tsunami heights along the Bengkulu coastal area. This indicates that an appropriate rectangular source for tsunami forecasting can be estimated from the W-phase inversion results and the scaling relationships (Blaser et al. 2010) with a rigidity of 3.0 × 1010 Nm−2, which is a typical rigidity for shallow thrust earthquakes.

4.3 Fault Model for the 2010 Mentawai Tsunami Earthquake (Mw 7.8)

The vertical components of the broadband seismic data from 13 stations are used for the W-phase inversion for the 2010 Mentawai tsunami earthquake (Fig. 5a). The solution of the W-phase inversion provides a seismic moment of 6.58 × 1020 Nm (Mw 7.8), thrust mechanism (strike = 328°, dip = 7.2°, and rake = 101°), centroid location (latitude 3.71° S and longitude 99.12° E), and centroid depth of 17 km. The vertical components of the observed and computed waveforms are compared in Fig. 5b. The fault length and width estimated from the scaling relationships (Blaser et al. 2010) are 103 km and 51 km, respectively (Fig. 5c). The estimated seismic moment is slightly smaller than the total seismic moment of 1.0 × 1021 Nm (Mw 7.9) estimated from tsunami waveform inversion by Satake et al. (2013). The centroid location placed the earthquake in the oceanic plate off the trench. However, the earthquake should have occurred as an underthrust event along the subducted plate interface (Satake et al. 2013). We discuss the accuracy of the centroid location in terms of the tsunami forecast in the discussion section.

Fig. 5
figure 5

Fault model estimated for the 2010 Mentawai tsunami earthquake: a station distributions (red triangles) used in the W-phase inversion, b concatenated observed (black trace) and synthetic (red trace) W-phase waveforms for vertical component, and c mechanism estimated by the W-phase inversion and the fault model estimated from the scaling relationships. The yellow star represents the centroid location of the 2010 Mentawai tsunami earthquake

Full size image

4.4 Tsunami Simulation for the 2010 Mentawai Tsunami Earthquake

The 2010 Mentawai earthquake generated a much larger tsunami than expected from the earthquake magnitude. Satake et al. (2013) surveyed the tsunami heights and inundation distances at seven regions in the Mentawai Islands, i.e. Lakau, Maonai, Asahan, Sabeu Gunggung, Muntei, Macaroni, and Tumalei, as shown in Fig. 6. Tsunami numerical simulations for the 2010 Mentawai tsunami earthquake were conducted using the estimated fault model with a length of 103 km and width of 51 km as shown in Fig. 5. Tsunami inundation simulations were conducted with the moving boundary condition in the finest grid areas shown in Fig. 7. The best computed results at seven regions, shown in Fig. 7 with red rectangles, are compared with survey data by Satake et al. (2013) in Fig. 8. The tsunamis were computed using six different slip amounts (or corresponding rigidities for the seismic moment of 6.58 × 1020 Nm) as follows: 12.25 m (μ = 1.0 × 1010 Nm−2), 8.75 m (μ = 1.4 × 1010 Nm−2), 8.17 m (μ = 1.5 × 1010 Nm−2), 7.66 m (μ = 1.6 × 1010 Nm−2), 6.81 m (μ = 1.8 × 1010 Nm−2), 6.13 m (μ = 2.0 × 1010 Nm−2).

Fig. 6
figure 6

Comparison of the averaged measured tsunami heights (black bar) by Satake et al. (2013) and averaged computed tsunami heights (red bar) using the fault model with the best slip amount of 8.17 m (corresponding rigidity of 1.5 × 1010 Nm−2) for the 2010 Mentawai tsunami earthquake

Full size image
Fig. 7
figure 7

Computational areas of tsunami inundation with the moving boundary condition (blue boxes) and plot areas in Fig. 8 at seven locations: Tumalei, Macaroni, Sabeu Gunggung, Muntei, Asahan, Maonai, and Lakau on Mentawai Islands (red boxes)

Full size image
Fig. 8
figure 8

Comparisons between the measured tsunami heights (red bars) by Satake et al. (2013) and those computed (blue bars) for the best model of the 2010 Mentawai tsunami earthquake at seven locations: a Macaroni, b Tumalei, c Sabeu Gunggung, d Muntei, e Asahan, f Maonai, and g Lakau. The computed tsunami inundation at each location is also shown

Full size image

The Aida numbers (K) were calculated for six slip/rigidity combinations mentioned above. For this calculation, the measured and computed tsunami heights at 23 measured points in seven locations (Fig. 8) were used. Figure 9 shows the calculated Aida number, K, with the slip amounts in the x axis. The best slip amount of 8.17 m (μ = 1.5 × 1010 Nm−2) was found for a K value of 1.00. The computed tsunami inundation areas and tsunami heights from the best slip amount are compared with the measured tsunami heights by Satake et al. (2013) as shown in Fig. 8 and Table 2.

Fig. 9
figure 9

Aida number K, the geometric mean ratio between the observed tsunami heights and those computed, with various slip amounts of the estimated fault model for the 2010 Mentawai tsunami earthquake

Full size image
Table 2 Comparison of the measured and computed tsunami heights for the 2010 Mentawai tsunami earthquake
Full size table

In the Macaroni area, the measured tsunami heights between 2.9 and 5.1 m are consistent with the computed heights of 3.7 to 4.1 m (Fig. 8a). In the Tumalei area, the measured tsunami heights between 4.6 and 5.4 m are slightly underestimated by the computed heights of 2.3 to 3.8 m (Fig. 8b). In the Sabeu Gunggung area, the measured tsunami heights from 5.0 to 5.8 m are similar to the computed heights of 6.1 to 6.6 m (Fig. 8c). In the Muntei area, the measured tsunami heights are between 4.6 and 8.8 m, and most are similar to the computed heights of between 4.8 and 6.3 m (Fig. 8d). However, the large measured tsunami heights of 8.8 m and 7.8 m are only slightly underestimated by the computed heights of 6.9 m and 5.9 m, respectively. In the Asahan area, the measured tsunami height of 7.4 m is slightly underestimated by the computed tsunami height of 5.7 m (Fig. 8e). In the Maonai area, the measured heights of 6.7 to 6.9 m are slightly overestimated by the computed heights of 8.0 to 8.4 m (Fig. 8f). In the Lakau area, the measured tsunami height of 2.5 m is consistent with the computed tsunami height of 2.0 m (Fig. 8g). These comparisons indicate that the measured tsunami heights are explained by the computed tsunami heights from the best model (slip = 8.17 m and μ = 1.5 × 1010 Nm−2). In Fig. 6, the average tsunami heights at the seven regions measured by Satake et al. (2013) are also compared with the computed heights. This also shows that the measured tsunami heights are explained by those computed from the best source model.

5 Discussion

The centroid location of the 2010 Mentawai tsunami earthquake using the W-phase inversion placed the earthquake in the oceanic plate off the trench (Fig. 5). Although we knew that the earthquake occurred along the subducted plate interface, the centroid location of the W-phase inversion off the trench was used to estimate the best fault model for the tsunami forecast to test how well the measured tsunami heights along the coast were explained. Because the technique presented in this paper is for tsunami forecasting of future large earthquakes, it is important to acknowledge that the result of the W-phase inversion may sometimes locate the earthquake seaward of the trench even when it actually occurs on the shallow megathrust. Our result shows that the tsunami computed from the best fault model, even though it is incorrectly located off the trench, agrees well with the measured tsunami heights along Mentawai Islands. This may indicate the robustness of our method in this study.

The 2007 great Bengkulu earthquake (Mw 8.4) was a typical subduction megathrust earthquake. The measured tsunami heights along the coast of Sumatra were well explained by the computed tsunami from the source model estimated by our method using the rigidity of 3.0 × 1010 Nm−2 (7.59-m slip amount). This rigidity is consistent with the depth-dependent rigidity curve estimated by Tanioka et al. (2017) at an estimated depth of 30 km (Fig. 10). The 2010 Mentawai tsunami earthquake (Mw 7.8) occurred at the plate interface near the trench. The measured tsunami heights along the coast of Mentawai Islands were also well explained by the computed tsunami from the source model estimated using the method with a rigidity of 1.5 × 1010 Nm−2 (8.17 m slip amount). This rigidity is slightly larger than the depth-dependent rigidity curve estimated by Tanioka et al. (2017) at an estimated depth of 17 km (Fig. 10). These results indicate that the rigidity differences along the plate interface as suggested by Bilek and Lay (1999) have a considerable impact on tsunami height estimation, as also suggested by Tanioka et al. (2017). The results also suggest that the blue dashed line in Fig. 10 should be used for the rigidity estimate along the plate interface off Sumatra in Indonesia. Therefore, appropriate slip amounts of source models for tsunami forecasts in this region can be rapidly estimated using the results of W-phase inversion and this depth-dependent rigidity (blue line in Fig. 10).

Fig. 10
figure 10

Comparison of the depth-dependent rigidity line (solid black line) suggested by Tanioka et al. (2017) and the rigidities estimated for the 2007 Bengkulu earthquake (green star) and the 2010 Mentawai earthquake (red star). The blue dashed line shows a depth-dependent rigidity line for the Central Sumatra region

Full size image

6 Conclusions

The fault parameters estimated from the W-phase inversion (Kanamori and Rivera 2008) and scaling relationships (Blaser et al. 2010) were used to obtain appropriate fault models for tsunami height estimates along the coasts resulting from the great 2007 Bengkulu earthquake and the 2010 Mentawai tsunami earthquake. For the 2007 Bengkulu earthquake, the measured tsunami heights by Borrero et al. (2009) were well explained by a single fault model using the rigidity of 3.0 × 1010 Nm−2 (7.59-m slip amount). For the 2010 Mentawai tsunami earthquake, the measured tsunami heights by Satake et al. (2013) were well explained by a single fault model using the rigidity of 1.5 × 1010 Nm−2 (8.17-m slip amount). These results show that it is important to use a depth-dependent rigidity relation adapted for local structures (such as the blue line in Fig. 10) in order to use the rectangular fault-based tsunami forecasting method of Tanioka et al. (2017) for Sumatra megathrust earthquakes.