A new dual earthquake and submarine landslide source model for the 28 September 2018 Palu (Sulawesi), Indonesia tsunami

ABSTRACT The September 2018 Palu (Sulawesi, Indonesia) tsunami has been a heavily debated event because multiple source models of three different types have been proposed for this tsunami: (i) The Mw 7.5 earthquake, (ii) landslides, and (iii) dual earthquake and landslide. Surprisingly, all of these three types of models were reported as being successful in the literature in terms of reproducing the existing tsunami observations. This can be partly attributed to the limited observations available for this tsunami. This study is motivated by the results of a marine bathymetric survey, which identified evidence for submarine landslides within the Palu Bay. Our modeling shows that the tsunami cannot be exclusively attributed to the Mw 7.5 earthquake. Inspired by the results of the marine survey, we propose a dual source model including a submarine landslide although most of the existing models include subaerial coastal landslides. Our dual model comprises an earthquake model, which has a length of 264 km, a width of 37 km, and a slip of 0–8.5 m, combined with a submarine landslide with a length of 1.0 km, a width of 2.0 km, and a thickness of 80.0 m located at 119.823°E and −0.792°S.


Introduction
The Palu city in Sulawesi of Eastern Indonesia was struck by a destructive tsunami following an M w 7.5 strike-slip earthquake on September 28, 2018 (Figure 1), which together killed more than 4,000 people. The earthquake origin time was reported as 10:02:45 (UTC) with a depth of 20.0 km and an epicenter at −0.256°S and 119.846°E by the United States Geological Survey (USGS, 2018). The Palu earthquake and tsunami produced the largest natural disaster in Indonesia following the 2004 Indian Ocean tsunami. Some other natural disasters in Indonesia in the aftermath of the 2004 event are the July 2006 Pangandaran earthquake (M w 7.7) and tsunami with over 800 deaths (Fujii and Satake 2006), the October 2010 Mentawai earthquake (M w 7.7) and tsunami with 408 fatalities (Satake et al. 2013, and the December 2018 Anak Krakatau volcano tsunami with over 450 casualties Mulia et al. 2020;Muhari et al. 2019).
Possibly, the 2018 Palu tsunami could be considered as one of the most debated tsunamis in the past few decades in terms of its source mechanism because different types of source models and in a large number have been proposed for this event, which sometimes contradict each other. Some of these models are given in Table 1. While some authors attributed the tsunami to only its earthquake source (an M w 7.5 strike-slip event), others completely ignored the earthquake source of the event and proposed that the tsunami was due to landslides (Table 1). Approximately four years after the 2018 Palu event and thanks to several field works and numerous modeling efforts, it is possibly fair to claim that there is a consensus in the tsunami community that the 2018 Palu tsunami was the result of a dual source, i.e. combined earthquake and landslide. There are several evidence for both co-seismic crustal deformation (e.g. field surveys by Natawidjaja et al. 2021) and submarine/subaerial landslides (e.g. marine surveys by Frederik et al. 2019; coastal surveys by Takagi et al. 2019;Aránguiz et al. 2020;Liu et al. 2020). Despite multiple efforts (Table 1) to characterize the earthquake and landslide sources of the tsunami separately, the topic is still open, and more research efforts are required as recommended by some authors listed in Table 1. Therefore, the purpose of this research is to offer a realistic dual-source model that is able to satisfactorily reproduce the existing tsunami observations (i.e. tide gauges and runup surveys) and address some of the shortcomings (discussed in the next Section) of the existing source models. We use data from the bathymetric survey of Frederik et al. (2019) to identify a potential submarine landslide source responsible for the Palu event combined with published seismic source models to propose our dual model. We apply numerical modeling to validate our source model through comparing tsunami modeling results with observed tide gauge data and surveyed runup heights.

A brief review of existing source models and our objectives
There have been extensive controversies in the past few years over the source mechanism of the 2018 Palu earthquake and tsunami, and multiple hypotheses have been proposed. Three different types of source models were proposed for this event (Table 1) comprising (i) Only-earthquake models (i.e. the coseismic crustal deformation), (ii) Only-landslide models (i.e. submarine or subaerial mass failures triggered by the earthquake), and (iii) Dual models (i.e. combined earthquake and landslide). Heidarzadeh et al. (2019) were among the first authors who showed that the tsunami was most likely the product of a combined earthquake and submarine landslide sources, and they approximated the location of a potential submarine landslide. Tsunami inversion by Gusman et al. (2019) revealed that an additional submarine landslide source is necessary to explain the observed inundation limits. Among the only-landslide models are those proposed by Pakoksung et al. (2019) and Nakata et al. (2020), who considered several hypothetical landslides and applied numerical modeling to show that some of them could be possible. On the other hand, Ulrich et al. (2019) and Jamelot et al. (2019) proposed models involving only earthquake ruptures. Among the authors who proposed dual sources are Heidarzadeh et al. (2019), Williamson et al. (2020) and Schambach et al. (2021) ( Table 1).
Several field surveys were conducted following the event to record coastal tsunami heights and damage Arikawa et al. 2018;Omira et al. 2019;Syamsidik et al. 2019). A marine bathymetric survey was conducted by Frederik et al. (2019)  following the event to map potential bathymetric changes in the Palu Bay. Takagi et al. (2019) and Liu et al. (2020) mapped several coastal landslides in the area through field surveys. All these field data of different types confirmed that the tsunami was most likely generated by a dual source, comprising the M w 7.5 strike-slip earthquake and landslides (submarine or subaerial). However, the current knowledge on a potential dual model for the 2018 Palu event is limited. There are three studies that proposed dual models: Heidarzadeh et al. (2019), Williamson et al. (2020) and Schambach et al. (2021). The study by Heidarzadeh et al. (2019) gives only the location of a potential submarine landslide without offering its dimensions or details. The two latter models are based on onshore subaerial landslides in their dual models with the exception that Williamson et al. (2020) also considered a hypothetical submarine landslide at the location previously proposed by Heidarzadeh et al. (2019) (Table 1). While there is evidence that such onshore subaerial landslides occurred (e.g. Takagi et al. 2019;Aránguiz et al. 2020;Liu et al. 2020), there is also evidence that submarine landslides were involved (e.g. marine bathymetric survey of Frederik et al. 2019;Heidarzadeh et al. 2019). Therefore, the shortcomings of the existing dual models are as follows: (i) Credible submarine landslides are missing in existing dual models although Heidarzadeh et al. (2019) proposed a hypothetical submarine landslide; and (ii) The simulated waves from landslide sources within the existing dual models show a mix of short and long waves as compared to the observed waveform in Pantoloan, and does not match the observations well.
Due to the relatively small sizes of coastal subaerial landslides, which result in some shorter-period waves compared to observations in Pantoloan, it is likely that a large submarine landslide was involved that could produce longer-period waves in Pantoloan. Marine surveys by Frederik et al. (2019) give evidence for such a large submarine landslide. The purpose of this research is to consider such a credible submarine landslide in our dual model. While we acknowledge that the dual-source models proposed by , Williamson et al. (2020), and Schambach et al. (2021) are important contributions, this study seeks to further complement them by considering a potentially reliable submarine landslide.
It is important to note that due to the limited tsunami observation data for the 2018 Palu event, many models of different types (earthquake, landslide, and dual models) can reproduce the observations. Therefore, the purpose of this study is not to discredit any of the published models. Rather, we aim at offering

Data and methods
The data used in this study comprise tide gauge records and runup heights of the tsunami, seismic source models for the 2018 Table 2). We acknowledge that multiple seismic source models were published for the 2018 Palu earthquake by different authors in the past few years (Table 1). Here, we use three of such source models in order to test a range of models and compare their performances in reproducing the observations. Table 2 gives the fault parameters for these three models, indicating that the maximum slip values for them are 8.5 m (USGS, 2018), 8.0 m (Jamelot et al. 2019), and 3.9 m (Wang et al. 2019). Here, the dislocation model of Okada (1985) is used to calculate the coseismic crustal deformation using the fault parameters listed in Table 2 for the three earthquake source models ( Figure 2a). The calculated crustal deformations ( Figure 2a) are used as initial conditions for tsunami modeling (e.g. Tsuji et al. 2011).
For defining the submarine landslide source models in this study, we use data from the actual post-event marine bathymetric surveys conducted by , it is very likely that the potential submarine landslide source was located around the middle of the bay because it is the location of the largest surveyed runup heights . Therefore, from the study of Frederik et al. (2019), here we only consider the potential landslide from the middle of the bay. In addition, we ignore the onshore subaerial landslides in our study by assuming that the effects of such onshore landslides are negligible as compared to those of a large submarine landslide. However, we acknowledge that multiple onshore subaerial landslides occurred during the Palu event (e.g. Takagi et al. 2019. Two transects of a potential landslide located around the middle of the Palu Bay (approximately 119.823°E and −0.792°S) are reported by Frederik et al. (2019) with a slope angle of 11-13° and thickness of 55-70 m. The two transects are distanced approximately 1,000 m. The length and width of the landslide cannot be precisely extracted from the report of Frederik et al. (2019) because the field data are limited; however, they can help us to approximate the length and width of the submarine landslide. Considering the bathymetric data of Frederik et al. (2019), we assume the length of the landslide in the range of 500-1,000 m and consider two scenarios LS-1 (length of 1,000 m) and LS-2 (length of 500 m) in order to consider the uncertainty of length estimation in our study and to conduct a sensitivity analysis (Table 3). In fact, a length in the range of 500-1,000 m appears meaningful given the results of bathymetric surveys of Frederik et al. (2019); however, the final length is determined by comparing simulated waveforms with observations. The width of the landslide is considered as 2,000 m for both LS-1 and LS-2 scenarios guided by the field data of Frederik et al. (2019). Inspired by the marine survey of Frederik et al. (2019), landslide thicknesses of 80 m and 40 m are considered for LS-1 and LS-2, respectively (Table 3).
To generate the initial sea surface displacement due to submarine landslides, we apply the empirical equations of Watts et al. (2005), which are based on studies previously conducted by Synolakis et al. (2002), Satake and Tanioka (2003), and Okal and Synolakis (2004). Figure 3a shows the initial sea surface displacement due to the two landslide scenarios, which are employed as initial conditions for tsunami simulations. Landslide tsunami modeling based on the empirical equations of Synolakis et al. (2002) and Watts et al. (2005) is a static modeling approach, which implies that the dynamic landslide generation process is overlooked. It has been shown by several authors that such an approach is successful in reproducing past landslide tsunami events (e.g. Synolakis et al. 2002;Tappin, Watts, and Grilli 2008;Satake 2015, 2017;Heidarzadeh et al. 2022b).
Tsunami modeling in this study is performed applying the numerical package COMCOT (Cornell Multi-grid   Comparison of observed and modeled tsunami waveforms at three tide gauge stations due to the three aforementioned earthquake source models. "U" and "S" denote "uplift" and "subsidence," respectively. "eps" indicates misfit 2 based on Equation (1). Misfits ("eps") are not calculated for Lahat as the tsunami signal is not clear at this station.  Comparison of observed and modeled tsunami waveforms at three tide gauge stations due to the two aforementioned landslide scenarios. "eps" indicates misfit 2 based on Equation (1). Misfits ("eps") are not calculated for Mamuju and Lahat as the simulated waves are very small at these stations. Coupled Tsunami model) (Wang and Liu 2006), which solves linear and nonlinear Shallow water Equations on Cartesian and Spherical bathymetric and topographic grids. Nonlinear equations were used in this study. We used the bathymetric and topographic grid provided by GEBCO (The General Bathymetric Chart of the Oceans), which has an original spatial resolution of 15 arc-sec (Weatherall et al. 2015). A single uniformly spaced bathymetric grid with grid spacing of 10 arcsec is used in this study, which is re-sampled from the original GEBCO grid in order to provide a smaller grid spacing for guaranteeing the stability of numerical computations. A grid size of 10 arcsec is equivalent to approximately 300 m (around earth latitude of zero). By knowing that the sea surface displacement field, according to the equations by Watts et al. (2005), is approximately three times of the landslide dimensions, the grid size of 10 arcsec (300 m) implies that there are at least nine or 10 grid points for each landslide source ( Figure 3a). Therefore, the resolution of the computational grid is sufficient to carry out the computations. A time step of 0.5 s was employed to guarantee the stability of the numerical scheme based on the Courant-Friedrichs-Lewy criterion (Courant et al. 1928). The total simulation time was 4.5 h. To avoid uncertainties associated with wave amplitudes during nonlinear wave inundation on dry land, in this study we used tsunami amplitudes at the 1-m water depth contour as approximations of tsunami runup heights (e.g. Tinti et al. 2006;Pranantyo et al. 2021).
Three types of simulations were conducted: (i) Simulations based on the earthquake source of the tsunami for three candidate earthquake models EQ-1, EQ-2, and EQ-3 as listed in Table 2 and shown in Figure 2a; (ii) Simulation using two landslide sources LS-1 and LS-2 (Table 3, Figure 3a); (iii) Simulations for dual (combined earthquake and landslide) sources, Dual-1 and Dual-2, where the earthquake and landslide sources are superimposed simultaneously (Table 4; Figure 4a). The basis for defining our dual models is explained in the next section. For modeling dual sources, it is assumed that the earthquake and the submarine landslide occur simultaneously.
To measure the qualities of fit between observations (O i ) and simulations (S i ), we use the misfit equation developed by Heidarzadeh et al. (2022a): where i ¼ 1; 2; 3; . . . :; N is a counter for the points on a time series, O i is an observation reading, S i is the corresponding simulation reading, and N is the total number of points in the time series. When simulations coincide with observations (a perfect match), 2 becomes zero; otherwise, it increases by an increase in mismatch between observations and simulations. We applied only the first tsunami wave recorded at each station for misfit calculations.

Modeling results and analysis
In this section, we compare modeling results with available observations of the 2018 Palu tsunami. Despite the great importance of the 2018 Palu tsunami, the existing observational data of the tsunami are limited. In the near-field, two tide gauge observations are available at Pantoloan and Mamuju ( Figure 1); however, the Mamuju record is outside of the Palu Bay and thus is most likely unaffected by the waves from the landslide source of the tsunami as the landslide waves are mostly confined within the Palu Bay ).
In addition, possibly the Mamuju record has a clock error of approximately 1 h. The Lahat Datu (abbreviated as "Lahat") tide gauge station, located approximately 600 km to the north of the epicenter, does not show any tsunami signal (Figure 2b). The other observation available for the Palu tsunami is data of tsunami runup height surveys conducted by several authors following the tsunami (e.g. Arikawa et al. 2018;Muhari et al. 2018;Omira et al. 2019). Before we discuss the results of simulations, it is helpful to note that, because the existing observational data of the Palu event are very limited, they can be reproduced by different types of source models. This is the reason that a large number of source models with different combinations have been proposed so far for this event (Table 1). Possibly, the truth about the source mechanism of the 2018 Palu event may not emerge until high-resolution bathymetric and seismic surveys, accompanied with seafloor coring of the Palu Bay, are conducted.

Performance of the earthquake source models (EQ-1, EQ-2, and EQ-3)
Comparisons of the simulated waveforms and coastal heights from the three earthquake source models are shown in Figures 2b and 5, respectively. Among the three earthquake models, Figure 2b reveals that EQ-1 performs better compared to EQ-2 and EQ-3 as the misfit resulting from EQ-1 at the Pantoloan station (eps = 0.28) is significantly smaller than that from EQ-2 (eps = 2.67) and EQ-3 (eps = 2.35). All the three Table 4. Information of two dual earthquake and landslide source models considered in this study for modeling the September 28, 2018 Palu (Sulawesi) tsunami.

Name of the dual model
Components of the dual model Earthquake model Landslide model Dual-1 EQ-1 LS-1 Dual-2 EQ-1 LS-2 earthquake models are unsuccessful in reproducing the Mamuju tide gauge record in terms of wave amplitudes and arrival times. Heidarzadeh et al. (2019) speculated that either the Mamuju tide gauge record has a clock error or the early tsunami waves were due to a secondary source such as a submarine landslide that occurred outside of the Palu Bay. Regarding the Lahat Datu station, all models produce very small waves. It is important to note that the EQ-1 model reproduces the initial phase and the period of the observation but fails in matching the amplitude. All three models significantly underestimate the surveyed runup ( Figure 5). These are strong evidence that indicates a secondary source (e.g. coseismic landslide) was most likely involved.

Performance of the landslide source models (LS-1 and LS-2)
Although it is clear that the 2018 Palu event was due to a dual source (combined earthquake and landslide), here we present the simulation results of candidate landslide models with the aim of understanding the nature of waves generated by submarine landslides. The simulation results for two candidate landslide scenarios (LS-1 and LS-2) are shown in Figures 3b and 6. The model LS-2 significantly underestimates the Pantoloan tide gauge record ( Figure 3b) and coastal heights ( Figure 6). The misfit resulting from LS-1 at the Pantoloan station is 0.24, whereas it is 0.39 for LS-2 ( Figure 3). However, the purpose of this section is not to rule out any of the LS-1 or LS-2 scenarios because although they may not reproduce the observations, it is likely that a combination of them with an earthquake model could reproduce the observations. This is discussed in the next section.

Performance of the dual-source models (Dual-1 and Dual-2)
For our dual model, we combine the earthquake model EQ-1 with the two landslide models LS-1 and LS-2 and produce two dual models: Dual-1 and Dual-2 ( Table 4). The basis for choosing EQ-1 for making the dual models is that EQ-1 is the only earthquake model that successfully reproduces the initial phases of the observed tsunami on the Pantoloan tide gauge (Figure 2b). Simulation results for the dual models reveal that the model Dual-1 slightly overestimates the observed tide gauge waveform ( Figure 4b) and slightly underestimates the surveyed runup heights ( Figure 6). However, its performance is far better than that of the model Dual-2. In terms of waveform match at the Pantoloan station, model Dual-1 gives a misfit of 0.45, whereas Dual-2 results in a misfit of 0.47 ( Figure 4). For runup, the misfits are 5.75 for Dual-1 and 7.68 for Dual-2 ( Figure 6). It is noted that our model Dual-1 underestimates some runup measurement points including a few locations at the southwest of the Palu Bay. By considering uncertainties associated with field measurements of runup heights and numerical simulations, we may conclude that model Dual-1 can approximately reproduce the observations.

Discussions
In addition to our final dual model, Dual-1 in Table 4 and Figure 4, three more dual models were proposed in the past few years by Heidarzadeh et al. (2019), Schambach et al. (2021), and Williamson et al. (2020) (Table 1). There are two differences between our model (Dual-1) and the other three modes: (1) our model is based on a credible submarine landslide (supported by an actual post-event bathymetric survey), while the other three models either only involve coastal subaerial landslides or include a hypothetical submarine landslide; (2) our model includes a single large submarine landslide, whereas the other models incorporate seven (Schambach et al. 2021) and 11 (Williamson et al. 2020) coastal landslides although the model by Williamson et al. (2020) also involves a hypothetical submarine landslide whose location was previously proposed by Heidarzadeh et al. (2019). Similar to this study, the dual models of Schambach et al. (2021) and Williamson et al. (2020) are also validated using tide gauge records and runup data. By comparing the results of modeling from our dual model, Dual-1 in Table 4 and Figure 4, with those of the aforesaid models, it is found that all dual models approximately equally reproduce the observational data, and it is hard to differentiate them in terms of quality of fit between simulations and observations. This may appear unexpected given the significant difference among these dual models. We speculate that such little differences among the results of different dual models could be attributed to the following: (i) Tsunami observations for the 2018 Palu event are sparse and scant and, thus, do not provide sufficient redundancy to distinguish among models; and (ii) The Palu Bay area is a very narrow and small basin (length and width of the Palu Bay are approximately 30 km and 7 km, respectively); therefore, different external forcings in a similar range (such as the three dual models discussed here) may not be distinguishable by this small basin.
From a natural hazard mitigation point of view, it is important to note that the largest tsunami heights and damage were due to the submarine landslide rather than the earthquake (Figures 5-6). For example, Figure 5 shows that the earthquake-only models produce maximum runup heights of approximately 2 m, whereas the observed runup heights were more than 10 m. In fact, the 2018 Palu catastrophe was mostly generated by the submarine landslide that was

Conclusions
We studied the source mechanism of the 2018 Palu (Sulawesi) tsunami through a numerical modeling approach accompanied with validation of the modeling results using observational data. The main findings are as follows: • Since limited observations are available for the 2018 Palu tsunami, they could be reproduced by multiple source models of different types (earthquake, landslide, and dual models). This is likely the reason that several source models are proposed and published for this tsunami. • The earthquake source of the tsunami underestimates the surveyed runup by an order of magnitude, which indicates that a secondary source (e.g. a landslide) was most likely involved. Therefore, we conclude that the tsunami was most likely the result of a dual earthquake and landslide source. • Although most of the existing landslide models for the Palu event are of the type of subaerial coastal landslides, we show that the existing tsunami observations can be approximately reproduced using a submarine landslide located inside the Palu Bay whose location is confirmed through marine bathymetric surveys.