Submarine landslide source modeling using the 3D slope stability analysis method for the 2018 Palu, Sulawesi, tsunami

Studies have indicated that submarine landslides played an important role in the 2018 Sulawesi tsunami event, damaging the coast of Palu Bay in addition to the earthquake source. Most of these studies relied on observed coastal subaerial landslides to reproduce tsunamis but could still not fully explain the observational data. Recently, several numerical models included hypothesized submarine landslides that were taken into account to obtain a better explanation of the event. In this study, for the first time, submarine landslides were simulated by applying a numerical model based on Hovland's 3D slope stability analysis for cohesive–frictional soils. To specify landslide volume and location, the model assumed an elliptical slip surface on a vertical slope of 27 m of mesh-divided terrain and evaluated the minimum safety factor in each mesh area based on the surveyed soil property data extracted from the literature. The soil data were assumed as seabed conditions. The landslide output was then substituted into a two-layer numerical model based on a shallow-water equation to simulate tsunami propagation. The tsunamis induced by the submarine landslide that were modeled in this study were combined with the other tsunami components, i.e., coseismal deformation and tsunamis induced by previous literature's observed subaerial coastal collapse, and validated with various post-event field observational data, including tsunami run-up heights and flow depths around the bay, the inundation area around Palu city, waveforms recorded by the Pantoloan tide gauge, and video-inferred waveforms. The model generated several submarine landslides, with lengths of 0.2–2.0 km throughout Palu Bay. The results confirmed the existence of submarine landslide sources in the southern part of the bay and showed agreement with the observed tsunami data, including run-ups and flow depths. Furthermore, the simulated landslides also reproduced the video-inferred waveforms in three out of six locations. Although these calculated submarine landslides still cannot fully explain some of the observed tsunami data, they emphasize the possible submarine landslide locations in southern Palu Bay that should be studied and surveyed in the future.

1 Introduction

Indonesia has recently experienced hardship because of earthquakes and tsunamis. Approximately 2.5 million people have been exposed to a tsunami with a 500-year return period, and a recent study suggested that Indonesia ranks third in the world for the proportion of the population exposed to tsunamis (Løvholt et al., 2014). A total of nine tsunamis have been observed since 2000, eight of which were associated with a significant earthquake, including the deadly Indian Ocean tsunami in 2004, according to the National Geophysical Data Center/World Data Service (NGDC/WDS, 2021). These types of natural disasters have caused many economic losses and fatalities in Indonesia. The country experienced limitations in reconstruction after the 2004 tsunami (Suppasri et al., 2015). Horspool et al. (2014) indicated that the annual probability of experiencing a tsunami similar to the aforementioned 3.0 m tsunami is ∼ 0.1 %–10 % across the entire Indonesian coast. The country also experienced another tsunami in December 2018 in Anak Krakatau, which killed 450 people (Muhari et al., 2019; Heidarzadeh et al., 2020).

The M_w 7.5 earthquake that occurred on 28 September 2018, near Central Sulawesi, Indonesia (0.256^∘ S, 119.846^∘ E), was caused by strike-slip faulting at shallow depths within the interior of the Molucca Sea microplate, which is part of the broader Sunda tectonic plate (USGS, 2021). The ruptured fault is a segment of a major active strike-slip fault zone named the Palu–Koro fault system (Liu et al., 2020). A major surface rupture with a 4 m maximum left-lateral offset was identified during a field survey in the NW Palu Valley (Patria et al., 2020). In addition, a 5.8 m surface rupture in Pawunu village occurred approximately 15 km south of Palu city (Badan Geologi, 2018). An unanticipated tsunami was triggered by this earthquake and caused damage around Palu Bay. The cumulative impacts from this tsunami led to approximately 2100 fatalities and more than 68 000 damaged houses (Pang et al., 2018; Hadi and Kurniawati, 2018). Several post-earthquake-tsunami field surveys reported damage along the Palu Bay coast and indicated run-up heights above 6.0 m (Mikami et al., 2019), specifically ∼ 9.1 m run-up heights and 8.7 m inundation depths at Benteng village (Omira et al., 2019) and ∼ 10.7 m run-up heights in Tondo (Widiyanto et al., 2019). Outside of the bay, there is little evidence of tsunami-flooding-related disruption, implying that most tsunamigenic sources were located within the bay (Liu et al., 2020)

1.1 Review of the attempts to model tsunami sources

Many attempts have been made to reproduce Palu tsunami events through modeling; however, these studies produced different outputs based on various tsunami sources. The finite fault model proposed by the USGS created only small vertical ground displacements that were not large enough to cause significant tsunamis within the bay (Pakoksung et al., 2019). Augmentations to the sources of the tsunami have been made. Due to the difficulty in reproducing the observed tsunami run-up heights and the Pantoloan tide gauge record (Heidarzadeh et al., 2019), landslide sources have been proposed as an additional source. At the Pantoloan tide gauge station in Palu Bay, a time series of a tsunami waveform was recorded. Several videos of the tsunami's arrival were also captured. Carvajal et al. (2019) used these videos to calculate the tsunami's arrival times and infer the waveforms at each site and emphasized that coastal/submarine landslides were the main contributing sources of this Palu tsunami.

Simulations of earthquake-generated tsunamis and submarine landslide-generated tsunamis involve different approaches. The generation of a tsunami induced by submarine landslides is a complicated process. Coupled dynamic models for earthquake-generated tsunamis have been well established and are commonly used, while coupled dynamic models for landslide-generated tsunamis are still limited (Heinrich et al., 2001; Pakoksung et al., 2019). Various modeling efforts have been carried out to simulate tsunami propagation using the COMCOT model (Heidarzadeh et al., 2019; Carvajal et al., 2019; Gusman et al., 2019; Liu et al., 2020; Sepúlveda et al., 2020). In more recent works, several studies have been performed to simulate landslide-induced tsunamis using other uncoupled methods, in which the landslide mass and tsunami propagation are separately calculated (Nakata et al., 2020; Paris et al., 2020; Ulrich et al., 2019).

Some previous works have successfully acquired agreeable results by comparing the surveyed run-up heights to the inundation depths, the recorded water levels at the Pantoloan tide gauge, or video-inferred waveforms. Many previous studies have used numerical models to reproduce a waveform at the Pantoloan tide gauge and have obtained good agreement with recorded data (e.g., Heidarzadeh et al., 2019; Carvajal et al., 2019; Takagi et al., 2019; Jamelot et al., 2019). Additionally, several studies have attempted to model the tsunami run-up heights and inundation depths, have compared them to the surveyed data, and have achieved acceptable comparisons (e.g., Pakoksung et al., 2019; Ulrich et al., 2019). Moreover, some studies have focused on inundation in Palu city. For example, Gusman et al. (2019) performed a simulation using vertical displacement and landslide sources to determine the inundation depths and area and compared the results with high-resolution field measurements by Paulik et al. (2019). More updated studies have attempted to reproduce video-footage-based waveforms, as proposed by Carvajal et al. (2019). In some works, hypothesized landslide sources are required to explain video-inferred waveforms (e.g., Nakata et al., 2020; Schambach et al., 2021). Moreover, combining seismic fault slip with observed landslide sources can also explain the wave characteristics based on video footage. For instance, Sepúlveda et al. (2020) combined landslide sources with coseismic ground displacement measured by interferometric synthetic aperture radar (InSAR) with different fault geometry configurations and archived the simulated tsunamis that were consistent with all the tsunami data, including surveyed run-up heights and video-inferred waveforms. These studies conducted simulations with different numerical models that used different sources of tsunamis or sources of landslide data to explain certain types of observational data. Table 1 summarizes the previous literature with respect to numerical models, tsunami sources, landslide data sources, and verification data.

Table 1 Summary of the main characteristics of the 2018 Palu event and tsunami reported in previous literature.

1.2 Review of submarine landslides in Palu Bay

Several studies have observed subaerial/coastal landslides or submarine landslides that can be visually observed through various methods. Small submarine landslide sources were identified from eyewitness accounts and field surveys (Arikawa et al., 2018). Sassa and Takagawa (2019) suggested that liquefaction induced coastal land collapse, which resulted in liquefied sediment flows and eventually led to a tsunami. Carvajal et al. (2019) used video footage and satellite images to identify visible coastal land collapses, and these landslide sources were used by later studies, including Gusman et al. (2019) and Takagi et al. (2019). Liu et al. (2020) studied coastal landslides by comparing pre-earthquake bathymetry to their surveyed post-earthquake bathymetry. These observed landslide locations are mapped and summarized in Fig. 1. Since actual submarine landslides are difficult to detect, hypothesized landslides have been proposed as an alternative in more recent works. Pakoksung et al. (2019) introduced four mega-sized landslides with volumes of ∼ 10–38 million cubic meters: three outside of Palu Bay and one in the southern Palu Bay part. Similarly, Nakata et al.'s (2020) study strongly emphasized a potential submarine landslide in the nearshore area near Palu city (Fig. 1).

Figure 1 Literature review of observed landslides: A–O represents the locations of observed coastal landslides, and the table summarizes landslides used/proposed in each study. The areas defined by yellow rectangles depict the potential submarine landslide areas proposed by Heidarzadeh et al. (2019), and the blue rectangle depicts the potential submarine landslide area proposed by Nakata et al. (2020).

1.3 Research objectives

Earthquake-induced tsunami modeling in previous studies cannot fully explain observational data. A better model for reproducing these tsunami impacts may be to reproduce tsunami events. Many past studies rely on hypothesized submarine landslides (Takagi et al., 2019; Pakoksung et al., 2019; Nakata et al., 2020) and suggest various potential submarine landslides within Palu Bay. This raises the question of whether exact landslides could actually occur in these areas and requires further investigation regarding these undiscovered submarine landslides. Therefore, more studies should focus on modeling submarine landslide sources. To the best of our knowledge, no existing studies have considered modeling submarine landslide sources based on slope stability analysis and the available observed soil data. The objectives of this study are the following.

1. Generate the potential submarine landslide using a sophisticated landslide model based on 3D slope stability analysis (which has never been performed according to the existing literature), also based on the existing observational soil data, and investigate whether the simulated submarine landslides match the observations or are located within potential areas suggested by past studies.

2. Reproduce all the field observations of tsunami records, including the run-up height, flow depth (or inundation depth) around Palu Bay, Pantoloan tide gauge waveform, video-inferred waveforms, and inundation depths and area in Palu city, with the parameters calculated by tsunami simulation based on developed landslide model.

2 Methodology

This section is composed of three parts. First, the landslide model is introduced; second, the numerical model used for tsunami propagation is introduced. Last, the data used for model settings and observational data used for model verifications are presented.

Figure 2 presents the workflow of this study regarding modeling procedures. Generally, there are two models: landslide and tsunami models. The former model considers the significant measured soil parameters as an input, including unit weight, cohesion, and internal friction angle. The modeling requires the initial water level and groundwater conditions according to the bathymetry. Similar to Nakata et al. (2020) and Paris et al. (2020), the landslide model assumes the initial shape of the landslide mass to be an ellipsoid. The geometry of landslides includes the maxima of the major and minor radii of the ellipsoid and the soil layer thickness. The seismic force is considered an additional sliding force in this study, which is discussed in a later section. Once the parameters are set, the model calculates the landslide mass (or landslide thickness) and the location of the landslide based on the minimum safety factor; these values are input into the tsunami model. The second part of the modeling procedure simulates tsunami propagation. This study adopts the same methodology as Pakoksung et al. (2019), details of which are explained in a later section. The tsunamis induced by the submarine landslide that were modeled in this study were combined with the other tsunami components, including the tsunamis simulated by USGS's finite fault model (Pakoksung et al., 2019) due to coseismal deformation generated by the Palu, Sulawesi, earthquake, as well as tsunamis induced by the observed subaerial coastal collapse depicted in Fig. 1. The reproductions of tsunami run-up height, inundation depth, Pantoloan tide gauge waveform, video-inferred waveforms at six locations, and inundation area in Palu city are determined to validate the observational data. The comparison index (Aida, 1978) is used to quantitatively compare the measured and calculated run-up heights and inundation depths. The input soil parameters and landslide geometry are repeatedly calibrated and adjusted until the comparison index reaches agreeable values.

Figure 2 Diagram of modeling procedures.

2.1 Landslide model

2.1.1 Model description

This study adopts a model based on Hovland's (1977) method as the landslide model. In Hovland's method, the collapsed soil mass cut off by the slip surface is decomposed into square or rectangular soil columns. Hovland's method defines the three-dimensional factor of safety as the ratio of the total available resistance along the failure surface to the total mobilized stress along with it (Ahmed et al., 2012). The safety factor formula for Hovland's method is shown in Eq. (1):

$$\begin{array}{}\text{(1)}& F={\displaystyle \frac{{\sum }_x{\sum }_y[cA+\mathit{\gamma }z\mathrm{\Delta }x\mathrm{\Delta }y\cos (\mathrm{DIP}\left)\tan \mathit{\varphi }\right]}{{\sum }_x{\sum }_y\mathit{\gamma }z\mathrm{\Delta }x\mathrm{\Delta }y\sin {\mathit{\alpha }}_{yz}}},\end{array}$$

where F is the safety factor, c is the adhesive strength, ϕ is the internal friction angle, A is the area of the slip surface, Δx and Δy are the soil mass length and width when the column is projected onto the xy plane, α_yz is the angle formed by the slip plane and the yz plane, DIP is the inclination direction of the slip surface, z is the depth of the soil mass, and γ is the unit weight of the soil mass. From the above, the formula can be simplified as the proportion of the resistance of the entire slipping soil mass, Q=ΣQ_i, over the sliding force, S=ΣS_i, as Eq. (2):

$$\begin{array}{}\text{(2)}& F={\displaystyle \frac{\sum Q_i}{\sum S_i}}.\end{array}$$

For stability over seismic activity, the design horizontal seismic intensity is introduced to consider the influence of the earthquake. The horizontal force due to the earthquake motion is assumed to be a part of the weight of the soil mass in the horizontal direction, as shown in Fig. 3a. The design horizontal seismic intensity is given as a vector k in the same horizontal direction as the gradient direction of the calculated element for each soil column. At this time, the resistance force (Q_i) and sliding force (S_i) of each soil mass are described as follows.

$$\begin{array}{}\text{(3)}& {\displaystyle }{Q}_i=[c{A}_i+\left\{n_i\cdot (e_g+k)(t_i\cdot u_i)(W-W_i)\right\}\cdot \tan \mathit{\varphi }]u_i,\text{(4)}& {\displaystyle }{S}_i=-[t_i\cdot (e_{\mathrm{g}}+k\left)\right(W_i\left)\right]t_i,\end{array}$$

where A_i is the area of the sliding surface, W_i is the weight of the earth column, and n_i and t_i are the unit vectors perpendicular to and parallel to the sliding surface, respectively. u_i is the direction in which the slip body slides and is determined at the same time as the slide body. e_g is the unit vector in the direction of gravity (Fig. 3b).

Figure 3 Diagram of soil mass components. (a) Design horizontal seismic intensity and direction of the soil mass. (b) The shape of the soil mass and directional vectors used for calculation.

Generally, the input soil parameters that are typically required for numerical modeling are soil unit weight, tensile strength, Young's modulus, Poisson's ratio, cohesion, and angle of friction (Kumar et al., 2016). In this study, the landslide model requires three essential soil parameters: the saturated condition of the cohesion, the internal friction angle, and the unit weight of the soil. In this study, landslide failure was analyzed by considering a safety factor > 1.

2.1.2 Measured soil data and landslide model settings

Observations of the submarine soil properties in Palu Bay do not exist in the literature. This research adapts soil data from field observations near the bay. The soil properties are defined by the pressure- and shear-wave velocity (V_p and V_s) profiles inferred from horizontal-to-vertical spectral ratio (HVSR) microtremor inversion over the soil depth, which was proposed by Cipta et al. (2020). Subsequently, V_s is translated to the standard penetration test (SPT) count (N) by the empirical equation proposed by Kirar et al. (2016) for all soil types: V_s= 99.5N^0.345. Overall, based on the measured soil properties, the study area is categorized into three main parts: northern, central, and southern zone (named as upper, middle, and lower Palu Bay respectively, as shown in Fig. 4). By the mechanism in the landslide model, the soil layers are categorized into two main layers: the sliding layer(s) and the base layer. Land failure is allowed on only the sliding layers (as shown by the stratum in Fig. 5). This research considers only two strata based on the observed soil properties in Fig. 4 and applies those from observation sites 504, 402, and 1001 for strata in the upper bay, middle bay, and lower bay, respectively. In addition, the soil properties are assumed to be saturated soil conditions, which significantly affect the result more than dry soil conditions in the landslide model.

Figure 4 Soil condition data for observation site nos. 1001, 402, and 504 and the Palu Bay zones delineated by the soil properties. The table includes the shear-wave velocity profile P (V_p) and S (V_s) parameters inferred from HVSR microtremor inversion, as well as the soil density over the soil depth and the soil thickness in each layer.

Figure 5 Diagram of the underground soil column structure in the model.

After determining the SPT N value, the saturated soil cohesion and internal friction angle are calculated from the empirical equations proposed by Kumar et al. (2016). The cohesion (C) can be computed from $C=-\mathrm{2.2049}+\mathrm{6.484}N$, for N ranges from 2–30 (for N < 2, we assume a random low C value). Moreover, the angle of internal friction (ϕ) is calculated from ϕ=7N for N < 4. ϕ= 27.12 + 0.2857N when N is in the range of 4–50. The values of these soil properties are tested to find good agreement with the observed tsunami run-up height and inundation depth and can be slightly different from the actual computed values from the above equations. Groundwater is assumed to exist, and the groundwater level is equal to the seabed elevation for all areas. Based on the observed water level at the Pantoloan tidal gauge, the initial water level is set to 2.3 m above mean sea level (MSL), representing high-tide conditions (Pakoksung et al., 2019).

Because of the unknown exact value of soil properties that can be estimated from only the stated empirical equations, this study considers a variation in soil parameters according to the available soil property data (Fig. 4). For example, the measured soil unit weight varies from ∼ 11.7–12.4 kN m⁻³ for stratum no. 1 (the first ∼ 3–8 m of depth from the ground surface) in lower Palu Bay and ∼ 13.2–15.0 kN m⁻³ for stratum no. 2. Regarding the landslide geometry, this study assumes an elliptical sliding surface. The geometry is composed of a major radius (the longer radius of the ellipse) and a minor radius. The major radius is always set equally to 2 times the minor radius. This study varies the major radius from 0.5 to 1 km. We perform ∼ 70 configurations based on these properties, and the final parameter settings for the landslide model are summarized in Table 2.

Table 2 Parameter settings for the landslide model.

2.1.3 Seismic force

To evaluate the impact of the earthquake, the design horizontal seismic intensity was proposed by the Eq. (5) (Noda et al., 1975):

$$\begin{array}{}\text{(5)}& \left|k\right|={\displaystyle \frac{\mathrm{1}}{\mathrm{3}}}{\left({\displaystyle \frac{a}{g}}\right)}^{\mathrm{1}/\mathrm{3}},\end{array}$$

where k is the design horizontal seismic intensity, a is the maximum value of the ground surface acceleration in centimeters per square second (cm s⁻²), and g is the gravitational acceleration, which is equal to 980 cm s⁻². At the Palu site, the ground surface acceleration is approximately 75 %–97 % of the gravitational acceleration (https://earthquake.usgs.gov/earthquakes/eventpage/us1000h3p4/map, last access: 23 January 2021), which means k= 0.28–0.33. In this study, a k of 0.30 is applied equally to all Palu Bay areas.

2.2 Tsunami model

This study adopted the same tsunami model as used in Pakoksung et al.'s (2019) study. A two-layer computational model, TUNAMI-N2, was developed to solve nonlinear shallow-water equations with two interfacing layers and kinematic and dynamic boundary conditions at the seafloor, interface, and water surface (Imamura and Imteaz, 1995). The mathematical model used in the TUNAMI-N2 code is made up of a stratified medium of two layers. The first layer is composed of a homogeneous inviscid fluid with a constant density ρ₁ that represents seawater, while the second layer is composed of a fluidized granular substance with a density ρ_s and a porosity μ that represents air. The calculated submarine landslides were considered in this layer. The mean density of the fluidized debris is assumed to be constant and equal to ${\mathit{\rho }}_{\mathrm{2}}=(\mathrm{1}-\mathit{\mu }){\mathit{\rho }}_{\mathit{\sigma }+}\mathit{\mu }{\mathit{\rho }}_{\mathrm{1}}$ (Macías et al., 2015). A more detailed explanation is given in Pakoksung et al.'s (2019) study. In addition, there is a useful review of landslide tsunami models provided by Heidarzadeh et al. (2014) for further study.

The motion of submarine landslide is modeled from the starting point – the original submarine landslide location from the landslide model – to the ending point, with the sliding plane inclined downslope. The modeling of landslide movement will end when the new bathymetric slope reaches the critical sliding slope which is 15.0^∘. All the submarine landslides are assumed to start moving at the same time, and the sliding movement is considered for 10.5 min after the earthquake.

2.3 Topographic and bathymetry data

Regarding the depth of the considered areas, this analysis used bathymetric data supported by Badan Informasi Geospasial (BIG), Indonesia. This database was developed prior to the Palu tsunami of 2018. BIG provided bathymetric and topographic data for the entire Palu Bay region as well as the adjacent continental areas. The domain was set to a 27 m bathymetric resolution in the sea, yielding a domain of 1155×810, which covered the entire Palu Bay area (Fig. 1), and a 1 m resolution for the nearshore Palu city with the size of 4.051×1.390 km², resulting in a 4051×1390 grid. A constant-grid tsunami simulation was solved at each time step of 0.01 s. It should be noted that the terrain data used in this study are adjusted by adding the constant of 2.3 m due to the lack of inundation area at high-tide conditions, as recommended by Pakoksung et al. (2019).

The USGS fault source with the dimension of 200 km × 30 km and a maximum slip of ∼ 9 m (https://earthquake.usgs.gov/earthquakes/eventpage/us1000h3p4/finite-fault, last access: 23 January 2021) was also used to model tsunami events. The initial water level was calculated by assuming that the alterations in the sea surface were equal to the seabed deformation based on a rectangular fault model

2.4 Comparison index

To quantitatively compare the measured and estimated run-up heights and inundation depths, Aida's (1978) correlation values; the geometric mean; and K and its variance, κ, are used and described by Eqs. (6) and (7) respectively:

$$\begin{array}{}\text{(6)}& {\displaystyle }\log K={\displaystyle \frac{\mathrm{1}}{n}}{\sum }_{i=\mathrm{1}}^n\log k_i,\text{(7)}& {\displaystyle }\log \mathit{\kappa }={\left[{\displaystyle \frac{\mathrm{1}}{n}}{\sum }_{i=\mathrm{1}}^n(\log k_i{)}^{\mathrm{2}}-(\log K{)}^{\mathrm{2}}\right]}^{\mathrm{0.5}},\end{array}$$

where $k_i=R_i/H_i$. R_i is the field measurement's run-up height at point i. H_i is the calculated run-up height at point i from the simulation, and n is the total amount of data. A K value close to 1 and small κ indicate good agreement between the observations and simulations.

3 Simulation results

3.1 Submarine landslides

This study configured ∼ 70 cases with variations in soil properties, angle of inclination, and earthquake intensity for the best estimation. A total of 23 submarine landslide events were simulated and are shown in Fig. 6 (landslide no. 1–23), excluding the proposed coastal/submarine landslides collected from past studies. The maximum landslide thickness varied from ∼ 5–20 m. Due to the stronger soil properties, fewer submarine landslides appeared in the northern bay, more land failure occurred in the central bay, and the southern bay zone was the most vulnerable to the submarine slope failure due to weaker input soil properties. The largest submarine landslide occurred around the nearshore area of Palu city, with a size of ∼ 0.04 km³, which is much smaller than the submarine landslides proposed by Pakoksung et al. (2019) and Nakata et al. (2020).

Figure 6 Simulated submarine landslides calculated based on Hovland's slope stability analysis in this study (landslide no. 1–23) and surveyed coastal/submarine landslides in previous literature (landslide A–O). The white and red color bar indicates the simulated landslide thickness for each grid (27 m resolution), and the black and white color bar indicates the observed landslide thickness. The dashed rectangles represent the potential submarine landslide areas proposed by Heidarzadeh et al. (2019), the dashed ellipse represents that proposed by Nakata et al. (2020), and the yellow ellipse represents that proposed by Schambach et al. (2021). The dots represent the video-inferred waveform observation points. The grey dashed line represents Palu–Koro fault.

However, most submarine landslides occurred in the potential areas proposed by past research. Heidarzadeh et al. (2019) proposed possible major submarine landslides along underwater slopes by using backward tsunami ray tracing, in which they placed a point tsunami source at the location of the Pantoloan tide gauge and propagated the tsunami model for 5 min (Fig. 6, dashed rectangle). Based on slope stability analysis, this study provided a maximum of eight submarine landslides that could occur at that location. The largest computed submarine landslide was located around the area suggested by Nakata et al. (2020), as shown in Fig. 6. The lower blue dashed ellipses hint at a large submarine landslide (Fig. 6, landslide no. 21) with a volume of 0.07 km³. Moreover, this study found three submarine mass failures with a total size of ∼ 0.05 km³ in the mentioned area. By following Nakata et al. (2020), Schambach et al. (2021) modeled submarine mass failure with a size of ∼ 0.026 km³ in a similar area and found good agreement with the nearby observed run-ups (Fig. 6, landslide no. 23 in yellow ellipse). According to these resulting landslides, this study strongly emphasizes the potential submarine landslides in lower Palu Bay, for example, landslide no. 21–23. In addition, the computed submarine landslides in the mentioned lower dashed ellipse area follow the possible erosion zone after the earthquake occurred when compared to the surveyed post-earthquake bathymetry based on Liu et al. (2020).

Unlike the study by Pakoksung et al. (2019), who suggested hypothesized submarine landslides in a northern area outside of Palu Bay, this study found only two small mass failures around these locations. This study did not generate any submarine landslide in the upper blue dashed ellipse zone around the Pantoloan tide gauge (Fig. 6), as proposed by Nakata et al. (2020). Moreover, the landslide model failed to simulate landslides in western Palu, as suggested by Carvajal et al. (2019) (zones I, J, K, L, and M) based on slope stability analysis and the available soil properties. However, Sepúlved et al. (2020) insisted on these contributions as the main sources. Therefore, it is necessary to consider the combination of these submarine landslides in previous literature and the present simulations to reproduce tsunami events.

3.2 Simulated tsunamis

3.2.1 Run-up heights and inundation depths

The resulting tsunamis in Figs. 7–11 were calculated based on the combination of submarine landslides from past studies (Fig. 6, landslide A–O) and this study's estimations (Fig. 6, landslide no. 13). Figure 7 shows the simulated tsunami run-up height and the comparison between simulated and observed values. This study compares the calculated run-ups with the observational run-ups from government reports and published research, including (1) the Geological Agency of Indonesia's (Badan Geologi, 2018) report, (2) the Agency of Meteorology, Climatology, and Geophysics (BMKG) of Indonesia, (3) Omira et al. (2019), and (4) Widiyanto et al. (2019). These studies' surveys cover all the study areas of this paper around Palu Bay. Overall, the simulated tsunami run-up heights vary from ∼ 3.5 to 13.5 m on the west coast and ∼ 4.1 to 11.0 m on the east coast. Large wave amplitudes are found at the shoreline of the western coast in middle Palu Bay and the eastern coast in lower Palu Bay, which indicate slight overestimation in these areas. The comparison index including the geometric, mean K, and geometric standard deviation κ is featured in Fig. 7b. According to the figure, overall, the simulation is overestimated as K= 1.18, which is > 1.0, and κ= 1.44 from 98 observations. A comparison with Widiyanto et al. (2019) gives the best results when K= 1.17 and κ= 1.39 from 27 observation points. The scatter plot in Fig. 7b indicates that most of the computed run-ups are in the range of ±2.1 m error. It should be noted that the measured run-ups from past research were adjusted to our terrain dataset; that is, 2.3 m was added (Pakoksung et al., 2019).

Figure 7 Simulated tsunami run-up heights resulting from the modeled submarine landslides based on Hovland's slope stability analysis and combined with past studies' observed submarine landslides. All the submarine landslides presented in the figure were used to model the tsunami run-up height. (a) Calculated maximum tsunami amplitude (center) with the maximum calculated tsunami amplitude (red line) at the shoreline of the western coast (left) and the eastern coast (right), as well as the measured run-up heights from Badan Geologi's report (blue triangles), BMKG (green circles), Omira et al. (2019) (crosses), and Widiyanto et al. (2019) (yellow squares). (b) Comparison between simulated and observed tsunami run-up heights by Badan Geologi's report (blue), BMKG (green), Omira et al. (2019) (red), and Widiyanto et al. (2019) (yellow), as well as the comparison index, K, and κ values. n stands for the total number of observations. The gray shading represents the error related to the over-/underestimation, i.e., ± standard deviation (∼ 2.1 m).

Figure 8 Simulated tsunami inundation depth (flow depth) with the modeled submarine landslides and observed submarine landslides from previous studies. All the submarine landslides presented in the figure were used to model the tsunami flow depth. (a) Calculated maximum flow depth (center) with the maximum calculated flow depth (blue line) at the shoreline of the western coast (left) and the eastern coast (right) and the measured inundated depths from Badan Geologi's report (red triangles), Omira et al. (2019) (crosses), and Widiyanto et al. (2019) (yellow squares). (b) Comparison between simulated and observed tsunami inundation depths by the Badan Geologi (2018) report (blue), Omira et al. (2019) (red), and Widiyanto et al. (2019) (yellow), as well as the comparison index, K, and κ values. n stands for the total number of observations. The gray shading represents the error related to the over-/underestimation, i.e., ± standard deviation (∼ 1.8 m).

The tsunami inundation based on simulated submarine landslides is also evaluated using the evidence of inundation surveyed by past studies, including Badan Geologi's report, Omira et al. (2019), and Widiyanto et al. (2019). The tsunami flow depth apparently follows a similar trend with the run-up heights but with more fluctuations. Generally, there are overestimated trends on both sides of the bay (Fig. 8a), and the simulation is slightly overestimated at K= 1.13, which is > 1.0. and κ= 2.15 from 62 observations according to Fig. 8b. The inundation depth results show large variances when compared to the surveyed inundations. This study calculates the inundation depths by the simulated run-up height subtracted by the topographic terrain with a 27 m resolution. The coarser resolution may inevitably affect the simulation results. This study also performs a tsunami inundation simulation on Palu city with a 1 m resolution terrain data, which shows a significantly better result. We discuss this result in the next section.

3.2.2 Inundation area and depths around Palu city

Palu city, where many buildings for residential, commercial, industrial, governmental, short-term stay, and other purposes are located, has been impacted by inundation ∼ 300 m from the coastlines. Paulik et al. (2019) conducted a high-resolution field survey to measure tsunami inundation heights at 371 building sites and indicated inundation depths ranging from 0.1 to 3.65 m, with a mean of 1.05 m and a standard deviation of 0.55 m. Moreover, the tsunami simulation results in this study can be slightly overestimated. The calculated inundation depths range from 0.12 to 3.83 m, with a mean of 1.30 m and a standard deviation of 0.50 m above the ground. A comparison between the simulated and observed tsunami inundation depths is shown in Fig. 9a. A total of 164 observation spots were selected in this study due to the limited availability of 1 m resolution terrain data. The comparison index indicates that K= 1.17 and κ=  1.93, which means that the depths obtained by the simulations slightly outnumber Paulik et al.'s (2019) observed inundation depths (Fig. 9b). Although the comparison was made at different scales, the high resolution of the data resulted in a higher accuracy of the inundation depths.

Figure 9 Tsunami inundation depths and the area around Palu city. (a) Simulated inundation depths and observation depths by Paulik et al. (2019) plotted over longitude. (b) Comparison between the simulated and observed tsunami inundation, as well as the comparison index, K, and κ values. n stands for the total number of observations. (c) Simulated tsunami inundation area (in color) and the observed inundation area retrieved from Gusman et al. (2019) (in grayscale).

The boundary of the tsunami inundation area was retrieved from Gusman et al. (2019), who interpreted the boundary from visible tsunami debris in satellite images. The comparison between the observed inundation area in Palu city is shown in Fig. 9c and indicates that the calculated inundation area is ∼ 1.195 km², while the selected observation of the inundation area is ∼ 0.747 km². This means that the simulation result is overestimated by ∼ 60 %.

Figure 10 Time series of sea level at the Pantoloan tide gauge with de-tided sea level data. The black line represents the simulated waveform based on submarine landslides and seismic faults, the blue line represents the simulated waveform based on seismic fault, and the red dot represents the sea level record at the Pantoloan tidal gauge.

3.2.3 Pantoloan tide gauge waveform

The record tsunami wave amplitude time series at the Pantoloan tidal gauge with de-tided sea level is depicted in Fig. 10. The first positive wave reached ∼ 0.20 m within the first 1–2 min and was followed by the first negative wave at ∼ −1.80 m within 5 min. The record tsunami wave peak of ∼ 1.95 m was reached at the tidal gauge within 6 min. Several past studies attempted to reproduce the waveform at the Pantoloan tide gauge. Although there was an argument that the reliability of the 1 min sampling interval of the tide gauge record did not sufficiently capture the shorter-period characteristics of tsunamis (Carvajal et al., 2019), many studies have shown agreeable results by reproducing the water level observations at the Pantoloan tide gauge. For example, the model of Heidarzadeh et al. (2019), using purely tectonic sources, performed backward tsunami ray tracing and modeled 5 min tsunami propagation to indicate two potential submarine landslide areas (Fig. 6) to match the waveform at the location of the Pantoloan tide gauge. Nakata et al. (2020) assumed a submarine landslide source near the tide gauge (Fig. 6, upper blue dashed ellipse) to explain the waveform characteristics. Moreover, Pakoksung et al. (2019) suggested that assumed large submarine landslides should be located north outside of the bay to generate the observed waveform characteristics at the Pantoloan tide gauge. Therefore, with calculated submarine landslides based on Hovland's slope stability analysis, this study supports the hypothesized area of Heidarzadeh et al. (2019) and shows that the landslides located in the southern part of Palu Bay affect the water level at the Pantoloan tide gauge. However, the results do not show the absolute accuracy of the peak tsunami wave. The proposed submarine landslide sources combined with the observed landslides generate a sooner arrival time and larger peaks of tsunami waves than the recorded data. The simulated waveform has an extreme negative peak of ∼ −5.0 m at 3 min and a positive peak of ∼ 3.0 m at 5 min after landslide occurrence (Fig. 10). The peak error percentage is used to indicate the performance of the simulation and can be determined by Eq. (8):

$$\begin{array}{}\text{(8)}& \mathit{\%}{\mathrm{Error}}_{\mathrm{peak}}={\displaystyle \frac{{\mathrm{Peak}}_{\mathrm{simulated}}-{\mathrm{Peak}}_{\mathrm{observed}}}{{\mathrm{Peak}}_{\mathrm{observed}}}}\times \mathrm{100},\end{array}$$

where Peak_simulated represents the peak calculated from the model simulation and Peak_observed is the peak of the tsunami wave from the observations. According to the equation, the simulation yields a peak error of 54 % and a 3 min shorter arrival time. It should be noted that this study tried to simulate only the observed submarine landslide sources (as was conducted in Lui et al., 2020), and the results did not show good agreement with the observed waveform.

Figure 11 Time series of sea level for the 2018 Palu tsunami at six locations: Dupa, the KN Hotel, Pantoloan, Talise, Wani, and West Palu. The black line represents the simulated waveform based on submarine landslides and seismic faults, the blue line represents the simulated waveform based on seismic fault, and the red line represents the video-inferred waveform retrieved from Carvajal et al. (2019). The simulated and observed waveforms are de-tided and adjusted to the same model's terrain datum. The location of each CCTV system is shown in Fig. 6.

3.2.4 Video-inferred waveforms

Past studies have shown some difficulties in reproducing CCTV-based tsunami characteristics based solely on coseismic sources for every available footage since they were introduced as an observational dataset (Carvajal et al., 2019; Sepúlveda et al., 2020). The observed tsunami waveforms inferred from video footage have short arrival times, and periods of ∼ 2 min resulted in short wavelengths (Carvajal et al., 2019). In general, by visual observation, the results of simulated video-inferred waveforms have even longer arrival times than the waveforms inferred from CCTV footage by judging the first peak of the tsunami wave. Three out of six video-inferred waveforms had comparable characteristics in southern Palu Bay, including in Talise, the KN Hotel area, and West Palu (Fig. 11), while the simulation results cannot fully explain the observations in northern Palu Bay, Pantoloan, and Wani. The simulation demonstrates that an overestimation of calculated landslide sources resulted in more devastating tsunami waves with higher peaks and longer wave periods in southern Palu Bay. The peak errors are 173 %, 23 %, −21 %, 60 %, −50 %, and 7 % for Dupa, the KN Hotel, Pantoloan, Talise, Wani, and West Palu, respectively. The time of the peak in the KN Hotel, Talise, Wani, and West Palu is approximately 0.90, 0.85, 1.00, and 0.21 min before the observations, while the peak time in the Pantoloan and Dupa occurs shortly after the observation, by 1.50 and 0.10 min. However, there is uncertainty regarding the timing of the video-inferred time series because the video footage was apparently captured promptly after the tsunami rather than earthquake shaking. The uncertainty may also result in ±15 to ±30 s of arrival time error in this study (Sepúlveda et al., 2020).

4 Discussions

Regarding the model parameters for simulating submarine landslide sources, we configured the model parameters, i.e., saturated soil unit weight, cohesion, internal friction of angle, groundwater condition, mean sea level, and earthquake intensity, as close as possible to the existing data. Slope failures can be found only in some steep slope areas inside Palu Bay. With the measured soil properties that were assumed as seabed property data, the model cannot comprehensively reproduce the observed landslides, which are located near Palu Bay, for both location and mass size. However, the study still can produce potential submarine landslides and reproduce tsunamis that are fairly comparable with the observation. Based on Hovland's slope stability analysis and measured soil properties, there is no slope failure in the mentioned nearshore or coastal areas. Despite modeling failure, this study can simulate submarine landslides in the potential areas proposed by Heidarzadeh et al. (2019), Nakata et al. (2020), and Schambach et al. (2021), which can imply that based on the supported theory, undiscovered submarine landslides in those areas are possible, especially a large mass failure with a size of at least 0.05 km³ in southern Palu Bay (Fig. 6, blue dashed ellipse). However, this study overproduced submarine landslide sources in southern Palu Bay according to tsunami simulation results. The run-up heights, inundation depths, waveforms at the Pantoloan tide gauge, video-footage-based waveforms, inundation area, and depths in Palu city all generally show overestimated trends. In particular, a comparison with observations shows an obvious overestimation of the simulation results of the southern Palu areas. These results imply that the real soil properties have a strong spatial variation. The soil properties in the simulation of southern Palu Bay, which was sampled near the Palu River, are probably weaker than the realistic conditions, resulting in many larger submarine mass failures than the observed failures. In addition, this study cannot reproduce the landslide sources around nearshore Palu Bay, where past studies strongly indicate that they mainly contribute to the Palu tsunami, i.e., zones I, J, K, L, and M (Carvajal et al., 2019; Sepúlveda et al., 2020) and zones B and C (Liu et al., 2020; Schambach et al., 2021), because these areas probably have stronger soil properties than those in the simulation. Notably, a landslide source in the nearshore area around the Pantoloan tide gauge, as proposed by Nakata et al. (2020) (Fig. 6, upper blue dashed ellipse), which still cannot be explained by this study, also needs to be addressed.

This study applied a constant earthquake intensity to the entire domain; however, in reality, earthquake intensity varies spatially. The peak intensity was located in southeastern Palu Bay according to the USGS. It should be noted that we tried to vary the earthquake intensity between the three study zones, with a maximum intensity of 76 % for upper Palu Bay and 80 % intensity for middle and lower Palu Bay, but the results did not show a significant difference from those when a constant intensity was applied.

There is controversy regarding the Hovland method. Hovland's slope stability ignores the resistant force acting on the surface and the friction between soil columns and considers only the force acting on the slip surface. This leads to a smaller safety factor when compared to the other methods (e.g., Janbu) (Ahmed et al., 2012); in other words, Hovland's method tends to give an overestimated output. However, this method is simple and suitable for calculation in a large area or a high-resolution terrain. There are alternative methods to tackle the problem of overestimating submarine landslides, such as the Bishop and Janbu methods. Hungr et al. (1989) stated that the simplified Bishop method offers a simple and efficient alternative that is applicable to a wide range of practical problems.

Further field observations of submarine landslide sources are required to confirm this study's results, especially for submarine landslides in the deep lower Palu Bay, as suggested by past studies and emphasized by this study. Field observations are also needed to confirm the soil property estimations used in the study area.

Further research should focus on the variation in the time of landslide initiation. This study assumed that all slope failures occur at the same time and start immediately after the earthquake, which may not reflect reality, as there is potentially a difference in the time at which a landslide starts to move. Adding coseismic sources or applying the results from the fault-slip model may help improve this study; Sepúlveda et al. (2020) showed that the combination of earthquake-triggered landslides and the inverted coseismic ground displacements measured with InSAR can reproduce the Palu tsunami with good agreement.

5 Conclusions

Rather than employing a trial-and-error approach to the generation of hypothesized landslide sources, this study applied a landslide model based on slope stability analysis for the first time to provide a better understanding of the potential submarine landslide-induced tsunami phenomenon. The simulated tsunami, which was triggered by modeled submarine landslide sources based on Hovland's 3D slope stability analysis with observational soil property data and configured landslide geometry, shows that the submarine landslide in southern Palu Bay is consistent with past studies; however, the magnitudes of simulated submarine landslides are overestimated, resulting in an overestimation of tsunami parameters, including run-up height and inundation depths, around Palu Bay and Palu city. This implies that the measured soil parameters inland of Palu Bay used in this study cannot fully explain the observed landslides. It is difficult to validate the simulation results according to all the available observational data, including tsunami run-up and inundation depth around Palu Bay and Palu city, the water level at Pantoloan tide gauge, and video-footage-based waveforms. The simulation results still do not fully explain the observational data, especially the flow depth across Palu Bay. Field observations of submarine landslides, especially in southern Palu, are needed to confirm the results of this study. In future studies, more accurate landslide simulations should be used by applying a more realistic configuration of the landslide thickness, soil parameters, and landslide geometry or dividing the study zones into zones with finer resolution. Additionally, the various landslide occurrence times should also be investigated, as strongly recommended by previous literature.