Submarine landslide source models consistent with multiple tsunami records of the 2018 Palu tsunami, Sulawesi, Indonesia

Previous studies have suggested submarine landslides as sources of the tsunami that damaged coastal areas of Palu Bay after the 2018 Sulawesi earthquake. Indeed, tsunami run-up heights as high as 10 m determined by field surveys cannot be explained by the earthquake source alone although the earthquake is definitely the primary cause of the tsunami. The quantitatively re-examined results using the earthquake fault models reported so far showed that none of them could fully explain the observed tsunami data: tsunami waveforms inferred from video footage and the field survey run-up tsunami height distribution. Here, we present probable tsunami source models including submarine landslides that are consistent with the observed tsunami data. We simulated tsunamis generated by submarine landslides using a simplified depth-averaged two-dimensional model. The estimated submarine landslide model consisted of two sources in the northern and southern parts of the bay, and it explained the observed tsunami data well. Their volumes were 0.02 and 0.07 km3. The radius of the major axis and the maximum thickness of the initial paraboloid masses and the maximum horizontal velocity of the masses were 0.8 km, 40 m and 21 m/s in the northern bay, and 2.0 km, 15 m and 19 m/s in the southern bay, respectively. The landslide source in the northern bay needed to start to move about 70 s after the earthquake to match the calculated and observed arrival times.


Introduction
The 28 September 2018 M w 7.5 Sulawesi earthquake that occurred 70 km north of Palu, Indonesia (US Geological Survey (USGS) 2018; Japan Meteorological Agency 2018), was a supershear strike-slip earthquake with a rupture velocity of 4.1 km/s (Bao et al. 2019). Syntheticaperture radar (SAR) analyses revealed that along-track horizontal displacement by the crustal deformation was about 5 m crossing from north to south in Palu Bay (Geospatial Information Authority of Japan 2018). A devastating tsunami was generated after this earthquake and produced damage around Palu Bay.
Field surveys of the impacted areas that began several days after the earthquake revealed that the tsunami locally reached heights of 10 m above sea level around Palu Bay Pribadi et al. 2018). Field surveys of coastal areas around Palu Bay (Omira et al. 2019;Widiyanto et al. 2019) revealed the overall picture of the tsunami height distribution in the bay. Figure 1a shows coastal tsunami flow depths and run-up heights reported by those surveys. Coastal collapses at multiple locations were also reported (International Tsunami Survey Team (ITST) 2018; Arikawa et al. 2018;Sassa and Takagawa 2019). Seafloor bathymetric surveys after the earthquake were also conducted (Frederik et al. 2019;Takagi et al. 2019).
A time-series waveform of the tsunami was recorded at the Pantoloan tide gauge station in Palu Bay ( Fig. 1b; Muhari et al. 2018). Numerous videos of the tsunami arrival were also recorded, and Carvajal et al. (2019) determined the arrival times and inferred the waveforms (wave periods and amplitudes) of the tsunami at each location based on those videos (Fig. 1c). Their results provided important clues to the source of the tsunami.
Various tsunami sources have been modeled based on seismic data [Heidarzadeh et al. (2019), using data from USGS (2018), Ulrich et al. (2019)], Sentinel-2 satellite optical data (Jamelot et al. 2019), Sentinel-1 satellite SAR images (Gusman et al. 2019), and hypothetical submarine landslides (Pakoksung et al. 2019). Although those studies compared calculated tsunami waveforms with the Pantoloan tide gauge record and tsunami heights determined by field surveys, a detailed comparison with tsunami heights over the entire bay has not been performed except for Ulrich et al. (2019), and little quantitative consideration has been given to the video footage. Therefore, a source model that can capture the entire observed tsunami data needs to be built.
In this study, we modeled the tsunami accounting for both fault models and submarine landslides to refine the tsunami source. We first simulated the tsunami using fault models (USGS 2018;Jamelot et al. 2019) and reassessed how well the calculated heights represent those determined by field surveys. We attributed differences  Pribadi et al. (2018); purple circles, Muhari et al. (2018)] and run-up heights [black triangles, Omira et al. (2019); orange triangles, Widiyanto et al. (2019)]. Waveform comparison points are locations where the tsunami waveform was inferred from the Pantaloan tide gauge (see b) or from video footage [see c; Carvajal et al. (2019)]. Landslide candidates (LSC1 to LSC6, dashed ellipses; see text) are also shown. The left and right graphs show flow depths and run-up heights on the western and eastern coasts, respectively, and black bars indicate average run-up heights over intervals of 0.002° latitude. b The observed tsunami waveform at the Pantoloan tide gauge: left, the two-day waveform including the earthquake origin time, and right, the waveform (with tidal contributions removed) during the 30 min following the earthquake. c Video-inferred tsunami waveforms from Carvajal et al. (2019). The waveform at the Pantaloan tide gauge (gray curve from b) is shown for comparison between the calculated and observed run-up heights to the inability of earthquake faulting alone to explain the tsunami source; we therefore extended our modeling to include submarine landslides, and explored sources that fit the video-inferred waveforms (Carvajal et al. 2019) and field survey heights. We prioritized the videoinferred waveform at Pantoloan in refining the source model because the video-inferred waveform appeared as an earlier arrival than the tide gauge waveform. Because the reported coastal collapses include only areas above the sea surface, it is unclear whether they continued to the seafloor, and if and by how much they contributed to the observed tsunami heights. Therefore, instead of setting multiple coastal and/or submarine landslide sources at once, we sequentially added landslide sources as necessary to explain the video-inferred waveforms. The purpose of this study is to reproduce the video-inferred waveforms and field survey heights as few tsunami wave sources as possible.

Numerical model of landslide mass and tsunami
A generation of tsunami by submarine landslide is a complex process, and various numerical models have been developed to analyze the propagation of the resultant tsunamis (see review by Heidarzadeh et al. 2014; Yavari-Ramshe and Ataie-Ashtiani 2016). Two practical models are known for modeling the landslide part of the landslide tsunami model: a viscous fluid model and a granular material model. The former viscous fluid models include Imamura and Imteaz (1995), Yalciner et al. (2014), and Baba et al. (2019). It has reproduced well the run-up heights of tsunamis generated by submarine landslides in application to the case of Papua New Guinea earthquake in 1998 (Imamura and Hashi 2003). Pakoksung et al. (2019) analyzed the 2018 Sulawesi earthquake using this type of model. The latter granular material models include Iverson and George (2014), Ma et al. (2015), Grilli et al. (2019) and Paris et al. (2020), whose models' landslide part governed by Coulomb friction is originated from Savage and Hutter (1989). Grilli et al. (2019) and Paris et al. (2020) applied this type of model to the case of Krakatau in 2018, demonstrating mass and tsunami behavior well. Grilli et al. (2019) also showed that tsunami waveforms at points far away from the source did not show much difference between the two models, dense Newtonian fluid model and granular material model.
In this study, we adopted the latter granular material model and the method similar to Paris et al. (2020). Titan2D (Pitman et al. 2003;Patra et al. 2005;Titan2D 2016) was used for the calculation of the granular mass material, and JAGURS (Baba et al. 2017) was used for the calculation of the tsunami propagation. Titan2D is a depth-averaged model that stably calculates the rheological motion of a continuum granular mass using real topographic data. JAGURS is a tsunami simulation model that can solve a depth-averaged nonlinear longwave equation.
The equations of the mass model used in this study are shown in (1)-(3), adding a buoyancy effect to the original Titan2D model: where t represents time, the x-axis and the y-axis are along the slope, and the z-axis is a direction perpendicular to the x-y plane. h is the mass thickness, v x and v y are the mass velocity, and g x , g y , and g z are the components of gravitational acceleration to each axis, respectively. φ int and φ bed are the internal and bed friction angles of the mass, and r x and r x are the curvature of the local basal surface. k ap is the coefficient that changes depending on the state of the active and the passive mass pressure, and is a function of φ int and φ bed . Note that k ap is negative if the mass is spreading. The description is omitted here [see the references of Savage and Hutter (1989), Pitman et al. (2003)]. The original Titan2D model is a dry avalanche model. We used the calculation code of Titan2D after multiplying the gravitational acceleration by the coefficient k b in the same way as Paris et al. (2020) to take into account the buoyancy of mass moving in water. Here, k b = (ρ s − ρ w )/ρ s . ρ s and ρ w are the density of soil in water and water. k b = 0.5 and ρ s ρ w = 2 in this study with reference to debris flow density in Denlinger and Iverson (2001).
The equations of the tsunami propagation used in this study are shown in (4)-(6), which are depth-averaged nonlinear long-wave equations. Although the tsunami calculations were performed in a latitude and longitude spherical coordinate system, equations in a Cartesian coordinates system are shown for simplicity: where t represents time, the x-axis and the y-axis are on a horizontal plane. η is the sea level, and M and N are flow rates in x and y directions. D is the total depth and D = h + η, h is the depth, g is the gravitational acceleration, and n is the Manning roughness coefficient. d is the input sea level displacement, which is set to equal to the sea bottom displacement in this study. The results of our mass motion calculations were given as seabed deformations of tsunami calculations. For stable calculation, the landslide mass motion is calculated in advance, as in Paris et al. (2020) and Baba et al. (2019), and it is used as input to the tsunami propagation calculation as the temporal deformation of the seabed. Figure 2 shows a conceptual diagram of the tsunami generation associated with mass motion. Here, seabed deformation was converted from universal transverse Mercator coordinates to latitude and longitude coordinates. The resistance between mass and water was ignored. This is because the details of landslide such as degree of interaction that is greatly affected by the shape of the mass and the degree of mixing with water have been uncertain in this event.

Conditions of landslide mass calculation
The numerical model of mass is discretized using the finite volume method and solved by the Godunov solver (Pitman et al. 2003). The initial shape of the flowable landslide mass is given as a paraboloid, and the mass motion is calculated according to its rheology (here, we selected a Coulomb-type rheology) under the driving force of gravity. The parameters to be given are limited to the radii of the major and minor axes (here, for convenience, r 1 and r 2 , respectively, in km), the orientation angle θ (counter-clockwise rotation from east), the maximum thickness h (m), and the horizontal location of the paraboloid, and the internal and bed friction angles of the mass, φ int and φ bed , respectively. φ int represents the frictional resistance of the mass to collapse; we used a fixed value of 30° because this parameter does not strongly affect the results (Ogburn and Calder 2017). φ bed represents the mobility of the mass, with smaller values representing higher mobility, i.e., faster downslope flow of the mass. The flowing mass tends to stop when the slope approaches this angle. Figure 3 shows the diagram of the calculation procedure of mass motion.

Conditions of tsunami calculation
The numerical model of tsunami propagation is discretized using the finite difference method and solved explicitly by the staggered-leaf-frog difference scheme (Baba et al. 2017). The spatial grid length for the calculation was set to 0.27 arc-seconds (about 10 m), the time interval to 0.04 s, and the prediction time to 30 min. Calculations were performed in the area between 0.52° S and 0.91° S and between 119.685° E and 119.91° E. The matrix size of the spatial domain for calculation was 3000 × 5200. Sea level fluctuations were assumed to reflect the seafloor deformations without delay at intervals of 10 s for 3 min after the mass started to move. The Manning roughness coefficient of friction was set to 0.025 in all regions uniformly. Landside boundary conditions were set to run-up. On land, the friction coefficient generally increases due to structures and coastal forests. Therefore, the result obtained may slightly be overestimated. Six output points were used to compare the calculated waveforms with the video-inferred waveforms and the waveform observed at the tide gauge station (Additional file 1: Figure S1). For comparison, we also used the JAGURS model to simulate the propagation of a tsunami generated by the earthquake alone. Vertical fluctuations due to horizontal motions on steep slopes were not quantitatively large (Heidarzadeh et al. 2019) and are omitted here. Because the local sea level was about 1 m higher than mean sea level at the time of the earthquake (Fig. 1b), we lowered the topography and bathymetry by 1 m when calculating the tsunami propagation, and restored the topography, bathymetry, and sea level after the simulation. Aida's (1978) correlation values for tsunami amplitudes: geometric mean, K, and its variability, κ , were used to quantitatively compare the observed and calculated run-up heights:

Comparison index
where k i = R i /H i , R i is the field survey height at the ith point, H i is the calculated height at the ith point, and n is the number of data. A strong match between observations and calculations is indicated by a small κ value and K close to 1. Calculated run-up heights were evaluated as the heights of the farthest west (on the western coast of the bay) or east (on the eastern coast) inundated points within every 0.002° of latitude, averaging the calculated mesh values.

Topographic and bathymetric data
We used DEMNAS topographic and BATNAS bathymetric data from Bandan Informasi Geospasial (BIG, the Geospatial Information Agency of Indonesia) for our calculations, because this dataset represents the bathymetry of Palu Bay more precisely than the global bathymetric data of the General Bathymetric Chart of the Ocean released in 2014 (GEBCO_2014; Weatherall et al. 2015;Heidarzadeh et al. 2019). We used 0.27 arc-second bathymetric data interpolated from the 6 arc-second Fig. 3 Diagram of calculation procedure of mass motion BATNAS dataset and 0.27 arc-second topographic data from the DEMNAS dataset. Both datasets were combined into a single land-sea topographic dataset (Fig. 1a) and cut in advance by paraboloid surfaces to form the soil landslide masses used in the landslide simulations. Artificial coastal structures in the Wani and Pantoloan ports were not considered; calculated waveforms in the simulations may therefore have slight errors related to this difference.

Video-inferred waveforms
Carvajal et al. (2019) inferred tsunami waveforms from video footage at six sites: Wani, Pantoloan, Dupa, Talise, KN Hotel, and West Palu (Fig. 1c). The arrival time of the tsunami is considered reliable at Wani, Pantoloan, KN Hotel, and Talise because uninterrupted video records of both the ground shaking due to the earthquake and the tsunami were obtained at these sites. At Dupa and West Palu, there may be a time lag of up to 2 min. Inferred amplitudes were minimum values and periods had errors of up to 2 s (Carvajal et al. 2019). The waveforms were determined at a tidal height of 1 m.

Tide gauge waveforms
Tide gauge data at the Pantoloan station in Palu Bay were also provided by BIG. For analysis, the tidal component was removed from the observations by subtracting the theoretical values calculated by superposition of trigonometric functions (e.g., Murakami 1981). The observed waveform is shown in Fig. 1b. About 5 min after the earthquake, a weak upward wave was followed by a large downward wave with an amplitude of at least 2 m and a wave period of 3-4 min (Heidarzadeh et al. 2019). Since the 1-min sampling interval of the tide gauge was long to sufficiently resolve the short period of the tsunami wave caused by a landslide, the peak of the tsunami wave at that location may have been even larger.

Field survey data
Run-up heights (Omira et al. 2019;Widiyanto et al. 2019) and flow depths (Pribadi et al. 2018;Muhari et al. 2018) were observed during field surveys (Fig. 1a). It is difficult to know whether the flow depth data includes splash heights. Therefore, we compared our calculated run-up heights with run-up heights (Omira et al. 2019;Widiyanto et al. 2019) alone. We used run-up heights including the tidal elevation at the time of the earthquake from the reports. Field survey data have multiple height values in close proximity at some locations. For comparison with calculated values, we averaged the observed run-up heights over intervals of 0.002° latitude (Fig. 1a).

Earthquake-only source models
We simulated tsunamis generated by the earthquake alone using the fault models of USGS (2018) Figure 4 shows our tsunami simulation results obtained with these two models. The results using both models show that the simulated waveforms at the Pantoloan station are similar to the results of Jamelot et al. (2019), and the period of the calculated waveforms is in good agreement with the observed tide gauge waveform. The amplitude of the calculated waveform by the Jamelot model is better agreement with the observed tide gauge waveform than that by the USGS model, considering the response of the tide gauge. The both models show good agreement with the arrival time of the video-inferred waveform at Wani and Pantoloan. However, the calculated run-up heights were low and not consistent with field survey heights, especially around the southern part of the bay. Even in the northern part of the bay, the results do not explain the video-inferred short-period waveform at Pantoloan. These results suggest that other sources such as submarine landslides contributed to the tsunami.

Submarine landslide source models
In this section, we consider submarine landslides as possible sources of the tsunami. We searched for hypothetical landslides that could represent the video-inferred waveforms. Because the tsunami arrival times are reliable at Wani, Pantoloan, KN Hotel, and Talise, we first assumed a tsunami source off Wani to reproduce the tsunami waveforms at Wani and Pantoloan, and a source off KN Hotel to reproduce those at KN Hotel and Talise. We then considered a tsunami source that would fit the waveforms observed at Dupa and West Palu, the sites with less reliable arrival times. Six locations were selected as submarine landslide candidates (LSC1 to LSC6, Fig. 1a).

Tsunami sources in northern Palu Bay
Here, we consider a landslide tsunami source that can reproduce the video-inferred waveforms at Wani and Pantoloan. At Wani port, a ship's crew reported that the polarity of the first wave was down, as the ship bottomed out shortly after the earthquake ( Assuming that a submarine landslide occurred immediately after the earthquake, and considering the tsunami propagation time, candidate landslide sources should be located to the north (LSC1 in Fig. 1a) or south of Wani (LSC3 in Fig. 1a). We tentatively set a wave source on the slope north of Wani at around 0.68° S (model A, LSC1 in Fig. 1a). Model parameters are reported in Table 1. The landslide mass descends to the west along the westward slope, pushing the sea surface up to the west and pulling it down to the east. This model can reproduce the rise of the water level observed in the video-inferred waveform at Wani, but under-predicts that at Pantoloan (Additional file 1: Figure  S2a). Calculated run-up heights are relatively high around 0.68° S, where the wave source is closer than at Wani, but do not match the field survey heights (Additional file 1: Figure S2b). We also tentatively set a wave source on the slope south of Wani at around 0.75° S (model B, LSC3 in Fig. 1a, Table 1). The result shows the calculated tsunami wave arrives at Pantoloan earlier than Wani, which is not consistent with the video-inferred waveform (Additional file 1: Figure S2a). Therefore, we set the wave source off the coast of Wani (LSC2 in Fig. 1a), assuming that the submarine landslide occurred 70 s after the earthquake (model C, Table 1). The volume was set to 0.02 km 3 ; r 1 0.8 km; r 2 0.4 km; h 40 m; and Φ bed 6. Figure 5 shows the results of model C. Figure 6 shows temporal change of spatial distribution and a cross-sectional time-series of the mass motion and tsunami. As the mass falls to the west, a large amplitude negative tsunami wave occurs on the shallow east side. After that, the wave spread in all directions. The waveform and run-up heights calculated using this submarine landslide source are consistent with the video-inferred waveforms and field survey heights in the northern part of the bay (Fig. 5). Here, the arrival time of the tsunami was determined by the position of the landslide, and

Table 1 Parameters of submarine landslide mass sources
The time that the landslide masses started moving, t 0 , is reported in seconds after the earthquake. All landslide masses were assumed to have Φ int = 30° (see text) a Reported latitudes and longitudes are for the centers of the landslide masses   Fig. 4 for source model C the period and amplitude were adjusted by the horizontal size, thickness and Φ bed of the landslide, respectively. This model predicts the maximum horizontal velocity of the front of the landslide mass to have been 21 m/s, as measured for a 5-m-thick location at the tip of the mass (arrows in Fig. 6b). The maximum gradient of the slope in Fig. 6b was 16°.

Tsunami sources in southern Palu Bay
Video footage recorded at Talise shows that a bore-like upward wave rushed to the shore about 2 min after the earthquake. Similarly, although the view of the sea was obstructed in the video recorded at KN Hotel, multiple cameras captured the tsunami striking the hotel about 2 min after the earthquake (Carvajal et al. 2019). Figure 7 and Additional file 1: Figure S3 show results for submarine landslide candidates LSC4, LSC5, and LSC6 (Fig. 1a) in the southern part of Palu Bay, referred to as models D, E, and F, respectively (Table 1). The volume was set to 0.07 km 3 ; r 1 2.0 km; r 2 1.5 km; h 15 m; and Φ bed 2 or 4. Of these three models, model E best fits the observations; the calculated waveforms match well the video-inferred waveforms at Talise and KN Hotel, and the calculated run-up heights are consistent with the field surveys. The arrival times calculated in model D and F are not consistent with those observed at Talise and KN Hotel. Here, the arrival time of the tsunami was determined by the position of the landslide, and the period and amplitude were adjusted by the horizontal size, thickness and Φ bed of the landslide, respectively. Although the calculated waveform at West Palu does not match the observed arrival time, the difference is within 2 min of tolerance. The characteristic of the observed waveform that the second wave is larger than the first wave is reproduced. Figure 8 shows temporal change of spatial distribution and a cross-sectional time-series of the mass motion and tsunami. The maximum horizontal velocity of the front of the mass is calculated to have been 19 m/s and the maximum gradient of the slope in Fig. 8b was 10° in model E.

Combined models
We considered two models that combine sources in the northern and southern parts of the bay. Model G combines the earthquake source (USGS model) and model E, whereas model H combines models C and E (Table 1); Additional file 1: Figure S4 and Fig. 9 show the results of both models: K = 1.13 and κ = 1.46 for model G, and K = 1.10 and κ = 1.42 for model H. Figure 10 shows the snapshots of the mass motion and the tsunami wave for model H. Although the calculated height is slightly smaller than the field survey height, model H reproduces the heights of the field surveys well with no bias between the eastern and western coasts. Although model G has the same accuracy as model H as a whole, the calculated heights on the east coast are small, especially near Wani and Pantoloan.

Tsunami simulation with dispersion effect
The dispersion effect was not considered in the tsunami simulations in previous sections. According to the index of Glimsdal et al. (2013): where τ is the index representing dispersibility, h 0 the water depth at the source, L the distance from the source, and λ the source size. τ = 0.05 when h 0 = 200 m, L = 2 km, λ = 1000 m as typical values of model H. Therefore, the dispersion effect is not significant in this case. The water depth at the tsunami wave source is shallow, and the distance from the source to the waveform comparison point is short whereas the tsunami associated with the landslide has a smaller source size than the earthquake fault, and generates short-period waves. In addition, since the purpose of this study is to focus on the waveform for several minutes after the arrival of the tsunami, it can be evaluated without considering the dispersion effect. In fact, the tsunami simulation results of model H with dispersion effect (using the calculation option of JAGURS) had little effect on the tsunami amplitude.   Aida's (1978) K and κ were calculated from these data and are reported for each coast as well as for all sites (black). The solid and dashed lines indicate the ratio of 1.0 between the calculated and observed run-up heights and the ratios of 2.0 and 0.5, respectively tsunami had passed, as observed by offshore buoys during the 2011 Tohoku-oki earthquake (Kawai et al. 2011). However, no such change as great as 1 m was observed (Fig. 1b). Therefore, the crustal deformations indicated by these models conflict with the observed tidal change. Although the USGS model has a smaller crustal deformation in Pantoloan than Ulrich et al. (2019) and the Jamelot model, it has at least about 50 cm. The calculated waveform at Wani by the USGS model only showed the sea surface level of 50 cm to 1 m 3 min after the earthquake, which is almost same as the level due to initial crustal deformation. As a result, the calculated tsunami amplitude at Wani by the model was small and the little runup on land was simulated. Therefore, the USGS model could not fully represent tsunami heights near Wani and Pantoloan whereas it had the reasonable heights of the entire bay by complementing submarine landslides in the southern Palu Bay as model G.
Although Sassa and Takagawa (2019) suggested that primary tsunami was coastal liquefied gravity flowinduced tsunami, the identification of submarine landslides and quantitative evaluation comparing with the observed tsunami data were not performed. By their report of the satellite images, no coastal collapse was reported in the area between Dupa and Talise, where field surveys reported run-up heights of about 10 m. We simulated a tsunami due to a coastal collapse near Dupa (model I, Table 1); the result shows the late arrival at Dupa and Talise and a limited area of tsunami runup near the source (Additional file 1: Figure S5) that is clearly different from the observed tsunami height distribution. Therefore, coastal collapses indicated by the satellite images alone cannot be the primary tsunami source explaining the video-inferred waveforms.
Regarding the potential source of the southern part of Palu Bay, Heidarzadeh et al. (2019) suggested the large submarine landslides along the submerged slopes in their discussion and identified the location of the potential submarine landslide based on backward tsunami ray tracing. Sunny et al. (2019) analyzed the video contents from the airplane and the boat in southwestern Palu Bay and estimated the unbroken wave-crest height of 8.2 m and splash height of 28.4 m. They also suggested the occurrence of the large-scale submarine landslides, comparing with the features of the tsunami waveform at laboratory experiments. Although the observed tsunami wave and the suggested submarine landslide in the southwestern do not explain the arrival time of video-inferred waveform in the eastern coast of Palu Bay from the horizontal location of sources, they may have a comparable impact to the tsunami heights in the western coast with that of our model.
Although a comparison of bathymetric data surveyed before and after the earthquake over the whole Palu Bay was reported, it was concluded that the resolution of the pre-earthquake data was not high enough to clearly identify the location of the submarine landslide (Frederik et al. 2019). The local survey in the southwestern part of Palu Bay was done, and the differences between bathymetric data before and after the earthquake were reported (Takagi et al. 2019). However, in our tsunami simulation based on the bathymetric change data (model J, Table 1), the result does not explain video-inferred waveforms and the distribution of tsunami run-up heights enough (Additional file 1: Figure S5). Here, the seabed deformation was input to the tsunami calculation over 30 s. Therefore, this wave source alone cannot explain the observed tsunami data.
The earthquake and coastal collapses (alone or together) cannot fully represent the observed tsunami data, and other wave sources are necessary; the submarine landslide sources based on the observed tsunami data are promising candidates as main tsunami sources.

Summary
We used a simplified depth-averaged two-dimensional model to simulate the 2018 Palu tsunami, and considered submarine landslide sources that could reproduce the first few moments of the tsunami as observed in video footage and the run-up height distribution determined by field surveys.
Although the earthquake is definitely the primary cause of the tsunami, the quantitatively re-examined tsunami simulation results by the earthquake fault models reported so far showed that none of them could fully explain the observational data. Especially, the earthquake source models in previous studies (Ulrich et al. 2019;Jamelot et al. 2019) seemed to have overestimate subsidence or uplift around the Pantoloan tide gauge station where considerable average tidal level change was not observed. The model by the USGS (2018) had small tsunami amplitudes at Wani and Pantoloan. Therefore, an even better model is needed to explain this tsunami with seismic faults.
On the other hand, hypothesizing a submarine landslide, a source model that can explain both the video waveform and the field survey height was determined. Of the various source models considered, a model with two landslides, one each in the northern and southern parts of the bay (model H), reproduces the observational data well. Estimated submarine landslide sources in the northern and southern parts of the bay had volumes of 0.02 and 0.07 km 3 . The radius of the major axis, the maximum thickness and the maximum horizontal velocity of the determined paraboloid landslide masses were 0.8 km, 40 m and 21 m/s in the northern bay, and 2.0 km, 15 m and 19 m/s in the southern bay, respectively. The landslide source in the northern bay needed to start to move about 70 s after the earthquake to match the calculated and observed arrival times.
The location of the landslide in the southern part of the bay suggested by Heidarzadeh et al. (2019) could be more specifically shown by matching with the observed tsunami data in this study. The present analysis showed the existence of a submarine landslide source consistent with observed tsunami data and the estimate of its size and location. As a future task, the impacts on the results when a more accurate seismic fault model is combined with a submarine landslide or when a submarine landslide tsunami generation modeling is improved using more accurate dynamic processes need to be investigated.
Additional file 1: Figure S1. Bathymetric (10-m contour interval) and topographic maps of the waveform comparison sites (red stars; see red squares in Fig. 1) for the tsunami calculations. Figure S2. Calculated tsunamis generated by landslide sources in the northern part of the bay. The results (models A, left; and B, right) are shown. a Comparison between calculated (red curves) and observed waveforms (black curves, video-inferred waveforms from Carvajal et al. 2019; gray curve, waveform at the Pantoloan tide gauge). b Distributions of the maximum calculated tsunami amplitude over the entire bay for each fault model. The color scale indicates the tsunami amplitude. Left and right graphs compare the calculated tsunami run-up heights (red curves) with those (black bars) determined by field surveys on the western and eastern coasts, respectively. Figure S3. Calculated tsunamis generated by landslide sources in the southern part of the bay. The results (models D, left; and F, right) are shown. Same as Figure S2 for source models D and F. Figure S4. Comparison of calculated and observed run-up heights. The heights on the eastern (blue) and western coasts (green) for model G are shown. Aida's (1978) K and κ were calculated from these data and are reported for each coast as well as for all sites (black). The solid and dashed lines indicate the ratio of 1.0 between the calculated and observed run-up heights and the ratios of 2.0 and 0.5, respectively. Figure S5. Calculated tsunamis generated by a coastal collapse and Takagi et al. (2019) (models N and O). Same as Figure S2 for source models N and O.