New High-Resolution Modeling of the 2018 Palu Tsunami, Based on Supershear Earthquake Mechanisms and Mapped Coastal Landslides, Supports a Dual Source

The Mw 7.5 earthquake that struck Central Sulawesi, Indonesia, on September 28, 2018, was rapidly followed by coastal landslides and destructive tsunami waves within Palu Bay. Here, we present new tsunami modeling that supports a dual source mechanism from the supershear strike-slip earthquake and coastal landslides. Up until now the tsunami mechanism: earthquake, coastal landslides, or a combination of both, has remained controversial, because published research has been inconclusive; with some studies explaining most observations from the earthquake and others the landslides. Major challenges are the numerous different earthquake source models used in tsunami modeling, and that landslide mechanisms have been hypothetical. Here, we simulate tsunami generation using three published earthquake models, alone and in combination with seven coastal landslides identified in earlier work and confirmed by field and bathymetric evidence which, from video evidence, produced significant waves. To generate and propagate the tsunamis, we use a combination of two wave models, the 3D non-hydrostatic model NHWAVE and the 2D Boussinesq model FUNWAVE-TVD. Both models are nonlinear and address the physics of wave frequency dispersion critical in modeling tsunamis from landslides, which here, in NHWAVE are modeled as granular material. Our combined, earthquake and coastal landslide, simulations recreate all observed tsunami runups, except those in the southeast of Palu Bay where they were most elevated (10.5 m), as well as observations made in video recordings and at the Pantoloan Port tide gauge located within Palu Bay. With regard to the timing of tsunami impact on the coast, results from the dual landslide/earthquake sources, particularly those using the supershear earthquake models are in good agreement with reconstructed time series at most locations. Our new work shows that an additional tsunami mechanism is also necessary to explain the elevated tsunami observations in the southeast of Palu Bay. Using partial information from bathymetric surveys in this area we show that an additional, submarine landslide here, when simulated with the other coastal slides, and the supershear earthquake mechanism better explains the observations. This supports the need for future marine geology work in this area.


INTRODUCTION
On September 28, 2018 at 6:02:45 PM local time (10:02:45 AM UTC), a Mw 7.5 magnitude earthquake struck Central Sulawesi, Indonesia, with the epicenter located approximately 70 km north of the city of Palu (USGS, 2018) ( Figure 1). The earthquake ruptured the Palu-Koro fault system, a predominantly strike-slip left-lateral fault (e.g., Socquet et al., 2019), along which large magnitude earthquakes have occurred in the past (Watkinson and Hall, 2017), two of which caused tsunamis in Palu Bay within the last century, in 1927 and 1968 (Prasetya et al., 2001). The 2018 rupture was supershear (Bao et al., 2019;Socquet et al., 2019;Fang et al., 2019) (i.e., it propagated faster along the fault than the local shear wave velocity), and the resulting ground motions caused widespread damage throughout the western Central Sulawesi region. Inland, the earthquake triggered landslides and induced considerable liquefaction that resulted in major destruction and numerous fatalities (Bradley et al., 2019;Watkinson and Hall, 2019;Miyajima et al., 2019). From eyewitness accounts and video evidence (Sassa and Takagawa, 2019), almost immediately after the earthquake, numerous coastal areas along Palu Bay experienced landslides (see locations of main ones marked LS-in Figure 1B), which were rapidly followed by destructive tsunami waves Carvajal et al., 2019). The earthquake, ground liquefaction, landslides, and tsunamis resulted in 4,340 fatalities and approximately 68,500 buildings were damaged or destroyed (BNPB, 2019).
Tsunami elevation time series were measured at two tide gauges, in the far-field at Mamuju (−2.66+ N, on W Sulawesi) and the near-field in Pantoloan Port (P in Figure 1B) (BIG, 2018). In the months following the event, international research teams conducted field surveys, recording earthquake and tsunami damage, landslides, and tsunami runup and inundation. Red dots in Figure 1B Widiyanto et al. (2019) (note only data labeled as runup in these references was used). Earthquake shaking, tsunami generation (particularly by coastal landslides; e.g., Figure 2), and various tsunami impacts were recorded in many amateur videos posted on social media, that were collected and analyzed (e.g., Carvajal et al., 2019), providing critical information on the timing and sequence of events. From the field and marine surveys, videos, and survivor accounts, it is clear that the earthquake, coastal landslides, and tsunamis closely followed each other, with major tsunami impact often taking place within minutes of the first shaking. Runups were highest in the south of the bay, reaching up to 10.5 m in the SE (Figure 3). At many locations where the coast was low lying, inundation only penetrated a short distance inland, which was interpreted as evidence of a landslide rather than an earthquake as the main tsunami generation mechanism (e.g., Muhari et al., 2018).
The only analogous recent events to 2018 Palu were Flores Island in 1992 and Gulf of Izmit in 1999. The Flores Island event was a shallow dipping thrust which triggered a coastal slide (Imamura et al., 1995), but there is no marine mapping to validate the slide size or volume. The Gulf of Izmit event is similar to Palu, with a strike-slip earthquake along the very active North Anatolian Fault, which triggered coastal landslides (Altinok et al., 2001). But again these landslides have not been mapped. This makes Palu an important event in that, for the first time, we have numerous tsunami data, including a comprehensive video data set, not only to fully investigate tsunami generation and coastal impact, but also to discriminate between the two, very different, source mechanisms, earthquake and landslide, and their contributions to tsunami hazard; and for the latter to confirm the importance of including dispersive effects in tsunami modeling. The improved understanding and modeling of the Palu event in this work can help identify, model, and more fully assess tsunami coastal hazards resulting from other similar tectonic environments.
Although there are now many tsunami simulations of the 2018 Palu event (e.g., Heidarzadeh et al., 2019;Takagi et al., 2019;Carvajal et al., 2019;Pakoksung et al., 2019;Gusman et al., 2019;Jamelot et al., 2019;Ulrich et al., 2019;Goda et al., 2019;Nakata et al., 2020;Sepúlveda et al., 2020;Liu et al., 2020, see summary of studies characteristics in Table 1), the tsunami mechanism, earthquake, coastal landslides, or both in combination, is still uncertain. In addition, whereas many published models simulate some, or even most, recorded runups around the Bay, the mechanisms are often ad hoc and do not reproduce the timing of tsunami waves from eyewitness accounts or the video evidence. Most coseismic tsunami sources (e.g., USGS, 2018;Socquet et al., 2019;Jamelot et al., 2019;Yolsal-Çevikbilen and Taymaz, 2019), are based on a primarily horizontal strike-slip earthquake mechanism, with limited vertical seabed motion (1-2 m). Theoretically, this should not be strongly tsunamigenic and should not, therefore, generate the elevated tsunami runups recorded from the south of the bay. There are many different interpretations of where the rupture is located under Palu Bay. In some earthquake models, the rupture crosses the bay as a simple, north-south, trending, connection (e.g., Socquet et al., 2019;Ulrich et al., 2019) (Figure 1B). In others, there is a change in direction under the bay (e.g., Jamelot et al., 2019) (Figure 1B), and some locate the rupture along the west coast (e.g., Song et al., 2019). When we completed our work, there was no resolution to these alternatives from the multibeam bathymetric data acquired by Frederik et al. (2019) in the deeper waters of the bay, because no seabed features had been identified as a possible rupture.
Hence, to quantify the effect on tsunami impact of this epistemic uncertainty in earthquake rupture, we simulated three representative coseismic sources: 1) Jamelot et al. (2019) and 2) Socquet et al. (2019) who inferred parameters for 9 and 294 sub-faults, respectively, by assimilating remotely-sensed observations of ground motion, and 3) Ulrich et al. (2019), who modeled the supershear seabed deformation as a function of space and time. Note that the latter more advanced study predicted a 1.5 m maximum vertical seabed motion.
Other coseismic mechanisms, derived from geodetic observations, yield larger vertical seabed motions and, therefore, could be more tsunamigenic (e.g., ∼ 3 m just south of the Balaesang Peninsula, Figure 1A, Song et al., 2019;Fang et al., 2019;He et al., 2019). These, however, are not within Palu Bay but farther north and, hence, cannot explain the tsunami here, particularly the fast arrival of large waves that impacted Palu City (the timing of events and waves will be detailed later).
Some have also argued (e.g., Ulrich et al., 2019) that the horizontal fault movement along the steep slope margins of Palu Bay resulted in an increased vertical water displacement causing elevated runups, in the manner proposed by Tanioka and Satake (1996). Hence, in our simulations of the three selected coseismic sources, we included this additional effect of enhanced vertical displacement as a function of the predicted horizontal fault movement.
The main challenge here, however, with the single earthquake mechanism, is that it cannot explain the timing of the tsunami impacts along the Bay from the coastal landslides reported in the survivor accounts and seen in video evidence, on land and in that captured by aircraft pilot Mafella flying over the bay shortly after the earthquake happened ( Figure 2).
With regard to the landslide tsunami mechanisms published so far for 2018 Palu (Table 1), Takagi et al. (2019) used a simplified numerical model of a dual earthquake/landslide source, with a single landslide located in the southwest of Palu Bay, mapped by high-resolution multibeam echosounder (MBES), whose tsunami was identified in the aircraft pilot video (Figure 2). Their model suggested that  Figure 1B), at tx108 s into the event (aircraft location in Figure 1B). "Boat" and "NBoat" mark where waves were also recorded on a small boat, as well as active subaerial slides (Supplementary Video S39 in Carvajal et al., 2019).  2019); (B) Footprint of BG with locations of (red dots) measured runups (Pribadi et al., 2018;Mikami et al., 2019;Omira et al., 2019;Putra et al., 2019;Widiyanto et al., 2019) (black dots) surface elevation time series inferred from shore-based videos (Carvajal et al., 2019, GM, grand mall; KN, KN Hotel; T, Talise; D, Dupa; P, Pantoloan; W, Wani) (yellow dots) observed landslides, (diamond) location of aircraft at 10:04:33 UTC, that filmed coastal landslides ( Figure 2). shorter period waves generated by the coastal landslide were followed by longer period waves from the earthquake. Pakoksung et al. (2019), Nakata et al. (2020), andSepúlveda et al. (2020) identified and modeled landslides as the most important, if not the principal, contributors to the tsunami. Their landslide parameters and corresponding wave generation, however, were hypothetical and selected to match observations. To date, only Liu et al. (2020) modeled landslide tsunamis from mapped landslide locations, but they: 1) did not model the tsunami generation from the landslides directly, using instead semi-empirical sources, and 2) did not simulate an additional earthquake mechanism. Finally all the landslide modeling studies to date simulated tsunami propagation with a non-dispersive model, which, as we will show, affected results.
Here, for the first time, we demonstrate that to explain the tsunami observations in Palu Bay requires simultaneously modeling both coseismic and landslide sources. We simulate the three coseismic sources discussed above (Jamelot et al., 2019;Socquet et al., 2019;Ulrich et al., 2019), the mapped (rather than hypothetical) landslides Takagi et al., 2019), and dual source combinations of these, using numerical models that include frequency dispersion effects Ma et al., 2012). We show that dispersion affects the shorter wavelength landslide tsunamis propagating into the deeper waters in the center of Palu Bay. We simulate the landslides as deforming granular material, with their tsunami generation, using the 3D physics-based numerical model of Ma et al. (2015). We use finer model grids and higher resolution bathymetric and topographic data in Palu Bay ( Figure 1) than in earlier work. We compare the results to a more comprehensive database of post-tsunami field survey results, including runups, tsunami elevations at the Pantoloan Port tide gauge together with those inferred at other locations from a novel analysis of the tsunami videos (Carvajal et al., 2019). From tsunami timing information in the aircraft pilot and other videos, we infer that there was a short delay in the triggering of the landslides by the earthquake, that we use in modeling and show that this improves the agreement with observations.
In Section 2, we detail and analyze tsunami observations, present the modeling methodology and data used to define tsunami sources and bathymetry/topography in model grids. In Section 3, we present model results for coseismic, landslide, and combined earthquake/ landslide tsunami simulations. Finally, in Section 4 we discuss the results and offer conclusions and perspectives for future work.

Tsunami Observations
In the following, we define t 0 as the start of the 2018 Palu event (10:02:45 AM UTC), i.e., the time the earthquake rupture begins at the epicenter (yellow star in Figure 1A).

Tide Gauge Data
Two operational tide gauges recorded apparent tsunami signals for the 2018 Palu event: 1) in Mamuju (−2.66+N, 118.89+E), in the Makassar Strait, on western Sulawesi about 250 km SSW from Palu Bay, a maximum trough-to-crest wave height of ∼ 0.25 m was recorded at t 19 min; and 2) in Pantoloan, within Palu Bay (−0.71157+N, 119.85731+E; site P in Figure 1B), a maximum trough-to-crest wave height of ∼ 3.8 m was recorded at t 5 − 6 min (BIG, 2018, Figure 4B). As noted in previous publications (e.g., Heidarzadeh et al., 2019), a tsunami wave traveling from the approximate location of the earthquake epicenter, north of Palu Bay ( Figure 1A) should take ∼ 45 min to arrive to the Mamuju tide gauge location, indicating either a clock error or that the signal here was caused by some other local source. In this work, as we focus on tsunami waves within Palu Bay, we do not use the  Okada (1985) and Tanioka and Satake (1996) NLSWE (Heinrich et al. (1998) and Hebert et al.
200 m resolution grid with two nested 10 m grids in palu city and pantoloan, derived from DEMNAS and BATNAS USGS (2018) and their own (hybrid source with nine subfaults parameterized from satellite data), using Okada (1985) and Tanioka and Satake (1996) for seafloor displacement None Ulrich et al.
Coupled EQ + tsunami model, seisol + StormFlash2D Triangular grid with maximum resolution 80 m in palu bay, derived from BATNAS Their own, modeled with seisol and coupled to the wave model, use Tanioka and Satake (1996)  Mamuju data nor try to explain this discrepancy. At Pantoloan, the pre-and post-tsunami tide gauge record shows that the earthquake did not cause measurable permanent changes in mean sea level (MSL) (Sepúlveda et al., 2020). Here, as noted by Carvajal et al. (2019), Sepúlveda et al. (2020), and Liu et al. (2020), the tide gauge is located in shallow water inside a harbor basin protected by a slotted seawall/dock (see Supplementary Figure S1 in supplementary material), which is not represented in the available bathymetric data nor in our model grids. While not affecting tide measurements, harbor structures may cause seiching and affect tsunami wave dynamics by modifying their elevation through reflection or dissipation; later in time ( > 650 s), the record may have been affected by waves reflected from the other side of the bay. The 1 Hz tide gauge measurements were averaged over 30 s and provided only every minute (BIG, 2018;Sepúlveda et al., 2020), yielding the fairly coarse time series plotted in Figure 4B. For this reason, while the gauge accurately measured long period waves, shorter waves, such as from landslides, were not recorded. This probably explains the difference in arrival time from the tide gauge data and closed circuit television (CCTV) video recording, overlooking the docks at Pantoloan Port (Supplementary Video S11 in Carvajal et al. (2019)), showing a train of shorter period waves arriving 2-3 min before the tide gauge registers any tsunami wave activity. In the following, the Pantoloan tide gauge data is compared to model results with these observations in mind.  Figure 1B. They discuss uncertainties in time series reconstructed from the CCTV videos and point out that while these are larger on surface elevation, due to the time stamp in the videos, phase information is quite good at all locations. This large video archive provided overwhelming evidence of tsunamis generated by coastal landslides (see https://agsweb.ucsd.edu/tsunami/2018-09-28_palu/carvajal_2019_videos_palu/). Most notably, a video of landslide tsunami generation on the western side of the bay was taken by Batik Airways pilot Ricoseta Mafella, at approximately 10:04:33 UTC, i.e.,t 108 s (Supplementary Video S38 in the archive, aircraft location marked by a magenta diamond in Figure 1B, at −0.829 + N, 119.869 + E; Mafella, personal communication), a composite picture of which is shown in Figure 2. The video shows multiple tsunamis generated as sets of concentric waves, offshore of the locations of coastal landslides LS-B, -C, -D, -E and -F* ( Figure 1B). Two of the smaller sets of waves at locations marked "Boat" and "NBoat" are consistent with a video made from a small boat at location "Boat" (Supplementary Video S39 in the archive); this video also shows the failure of a subaerial coastal landslide. Furthermore, a video taken from a ship docked at Taipa, on the southeast coast of the bay (marked in Figure 1B; Supplementary Video S31 in the archive) showed at least one other landslide failure to the north (potentially at location marked by LS-L or -M in Figure 1B). Sunny et al. (2019) analyzed the waves generated on the western side of the bay at locations LS-D, -E, and -F* ( Figure 1B), seen in the pilot Supplementary Video S38, using Google Earth to match camera viewing angles, and compared the observed wave measurements to dimensions of known objects. They estimated the widths of the sharp crescent waves to be 343 and 461 m, and the elevations (trough-to-crest) to be 24.1 m and 28.9 m at locations LS-D and -E, respectively ( Figure 1B). Based on the boat Supplementary Video S39, they estimated that at location -F*, the splash of the outgoing wave was 28.4 m high, and the elevation of the unbroken outgoing wave traveling toward the Boat location was 8.2 m. These estimated wave heights are all much greater than the initial values, on the order of 2-10 m, predicted by Liu et al. (2020) using semiempirical methods for the landslide tsunami waves generated at these locations.

Timing and Wave Sequence Analyses Based on Videos and Supershear Velocities
The combination of video evidence and supershear earthquake travel time can be used to estimate the time after rupture initiation that ground shaking began at various locations around Palu Bay.
CCTV Based on the consistent travel times and modeling estimates for the beginning of ground motion, we conclude that all locations within the bay most likely started shaking at t 9 − 15.5 s, which allows to constrain the tsunami wave arrival times from the videos that also show the start of earthquake shaking. At Wani and Pantoloan, therefore, we included a 9 s delay, to allow for seismic waves to reach this area, and in the southern sites of Dupa, Talise, KN Hotel, and Grand Mall, a 14 s delay. As mentioned above, at 10:04:33 UTC, or t 108 s, Pilot Mafella and his aircraft were located at −0.829 + N, 119.869 + E (near the eastern side of the Bay; Figure Figure 4A shows the short time series of surface elevation they estimated at the Wani dock (−0.6933+N, 119.8418+E), with an elevation wave cresting at 2 m. Crew members on the Sabuk Nusantara vessel, docked at Wani, reported that immediately after the shaking there was a ∼ 7 m withdrawal of the water, followed 3-5 min (180-300 s) later by a ∼ 15 m height wave cresting at ∼ 8 m (VOA-News, 2018). Although this estimated crest elevation is larger than based on the videos at the house, its timing is consistent with that of Carvajal et al.'s; additionally, the ship observation confirms there was a large depression wave (trough) preceding the crest.
• Supplementary Video S11, taken in nearby Pantoloan Port by a CCTV camera looking toward a crane on the dock (−0.7106+N, 119.8552+E), captured the initial tsunami waves at this location. Assuming that the video footage begins at the time of shaking, the trough of the initial shoreline withdrawal occurs at t 189 s, followed by a large positive wave at t 215 s. Figure 4B shows the time series of surface elevation estimated at this location (solid black line), with a −2 m trough followed by a 2.5 m elevation wave. This is consistent with waves inferred from the video recorded in Wani, but at the Pantoloan tide gauge, those shorter and higher waves were filtered out by the gauge (dashed black line).
• In Supplementary Videos S29-S31, taken in Taipa (−0.7794+N, 119.8580+E; Figure 1B), the timing of the videos is unknown and the time series of surface elevation could not be estimated. However, in Supplementary Video S29, Carvajal et al. (2019) note that there is a wave to the north that appears similar to other landslide generated waves located in other parts of the bay, which could potentially be attributed to sites LS-L or M ( Figure 1B).
• Supplementary Video S14, at Dupa (−0.8204+N, 119.8811+E; D in Figure 1B) begins some time after the earthquake shaking. Carvajal et al. (2019) estimated that the sea withdrawal started at t 105 s and Figure 4C shows the short time series of surface elevation estimated here, with a ∼ −1.5 m trough.
• In Supplementary Video S13, at Talise (−0.8589+N, 119.8789+E; T in Figure 1B), the water begins to withdraw at t 39 s, followed by a large wave striking the shore at tx39 s, as confirmed by the people transitioning from walking to running away from the coast. Figure 4D shows the short time series of surface elevation estimated here, with a −1.3 m trough followed by a ∼ 2 m crest.
• Six CCTV cameras were operated at the KN Hotel (−0.8650+N, 119.8775+E; KN in Figure 1B), ∼ 750 m south of Talise. The camera timestamps were adjusted to the time shaking started, at t 106 s based on Ulrich et al. (2019). A sea withdrawal is not seen, but tsunami inundation from a northerly direction begins at t 106 s. In Carvajal et al. (2019)'s Supplementary Video S3, the camera angle from the KN Hotel points across the bay in the direction of the LS-F* landslide. In Figure 5, showing video frames, a disturbance becomes visible at t 52 s (video time 38 s) in the upper right corner behind a tree. This could potentially be the wave generated by the LS-F* landslide. Figure 4E shows the short time series of surface elevation estimated here, with a ∼ 2 m crest.
• Finally, many videos were made from various floors of a parking structure in Palu Grand Mall (−0.8836+N, 119.8437+E; GM in Figure 1B), which were combined into a 11′20″ video  Figure 4F shows the short time series of surface elevation estimated here, with two waves with a −2.0 and 2.5 m largest trough and crest, respectively. dat; urls in their paper wrongly included a period or space after the word " Fig"). They estimated an approximate landslide volume of V S 3.2 10 m. Frederik et al. (2019) surveyed areas deeper than 200 m within Palu Bay, as well as outside of it, southwest of the Balaesang Peninsula. Within the bay, they could not find any clear sign in the bathymetry of a fault trace at seabed, nor any evidence that large deeper water submarine mass failures (SMF) occurred. Liu et al. (2020) surveyed the shallow coastal waters of Palu Bay and identified 14 locations where recent coastal landslides occurred. They estimated slide parameters based on differences between their new and the Badan Informasi Geospasial's (BIG; Geospatial Information Agency, Indonesia) pre-earthquake bathymetric contours. In a study published after our work was completed, Natawidjaja et al. (2020) reanalyzed the published bathymetries (Frederik et al., 2019;Liu et al., 2020) and interpreted the margins of the deep central channel in Palu Bay to be faults that were activated during the 2018 earthquake. They proposed several meters of vertical displacement for these, although this movement is within the m resolution of the data. These authors also suggested that the faulting triggered "massive" landslides in the southeast of the Bay, although comparisons of pre-and post-event data in this region by Liu et al. (2020) suggest that only the southern landslide was reactivated at this time.

Post-tsunami Surveys
In this work, we studied and modeled a subset of the coastal landslides clearly identified by Liu et al. (2020), labeled LS-B,-C,-D,-E,-F*,-L, and -M in Figure 1B ( Table 2), for which video evidence confirmed that wave generation occurred, using dimensions and volumes adapted from Liu et al. (2020). For this reason, with the exception of landslide F*, for which we used data provided by Takagi et al. (2019), we used the same labeling scheme as in Liu et al. (2020).
Tsunami Coastal Impact: Post-tsunami surveys were conducted by various international teams, in which flow depth and runup ( Figure 1B) were measured around Palu Bay.

Materials and Methods
We numerically simulate the 2018 Palu tsunami generation and propagation from earthquake or landslide sources, and the two in combination. We follow a methodology similar to that used in recent work by the authors and collaborators for other dual earthquake/landslide mechanisms (e.g., Tappin et al., 2014;Grilli et al., 2015;Grilli et al., 2017;Grilli et al., 2019;Schambach et al., 2019;Schambach et al., 2020). Using the best available bathymetric/topographic data, together with higher resolution computational grids than used in previous studies, we apply two state-of-the-art dispersive wave models: 1) the 3D non-hydrostatic wave model NHWAVE (Ma et al., 2012), with an underlying slide layer (Ma et al., 2015;Kirby et al., 2016), to simulate both the landslide motion and related tsunami generation, and 2) the 2D fully nonlinear Boussinesq wave model FUNWAVE-TVD  to simulate the propagation to the far-field and the coast of the superposition of landslide and coseismic tsunamis (or each individually), in nested grids of increasing resolution toward the shore. The grid data and modeling methodology are detailed next.

Label
Lon The study area encompasses Palu Bay (Figure 1), which is approximately 30 km long by 7.25 km wide with depths reaching up to ∼ 830 m ( Figure 6). Analyses of pre-and post-earthquake satellite images show either N-S or NNE-SSW strike-slip ground motion north of Wani, and NNW-SSE motion south of the Palu Grand Mall (e.g., Valkaniotis et al., 2018;Socquet et al., 2019), indicating that the rupture trace changed direction somewhere in the bay. However, as mentioned above, the initial interpretations of the high-resolution deeper water bathymetric data of Frederik et al. (2019), and the separate high resolution bathymetric survey by Liu et al. (2020) did not show any evidence of a clear rupture trace in the bay. In their recent interpretation of this data, Natawidjaja et al. (2020), in contrast, suggest several meters of vertical fault movement. Notwithstanding these new interpretations and, as noted by Liu et al. (2020), since the Palu-Koro fault was not previously mapped underwater, the existing earthquake models have adopted various assumptions regarding where the fault trace is located. Figure 1  Three Cartesian computational grids are used in tsunami simulations ( Figure 6; listed grid center coordinates are used as Mercator transverse geographic projection origin). The 30 m resolution base grid BG covers Palu Bay ( Figure 1) and is used in FUNWAVE and NHWAVE simulations; the NHWAVE BG grid also includes five boundary fitted, equally-spaced, layers in the vertical direction (from ocean surface to seabed). Two 7.5 m resolution grids are nested within BG, and used in FUNWAVE to more accurately model tsunami coastal impact in two areas of particular interest: 1) in the south of the bay, south grid SG includes the observation points of Grand Mall, KN Hotel, Taipa and Dupa; and 2) near and around Pantoloan, east grid EG includes the observation points of Wani and Pantoloan ( Figure 1B).
A variety of bathymetric and topographic data sets have been used in earlier modeling studies of this event. Here we similarly combine and interpolate onto our grids the best available bathymetric and topographic data to date for our study area. The resulting bathymetry and topography are shown in Figure 6. As an overall coarser data set, the earlier study of Heidarzadeh et al. (2019) used the General Bathymetric Chart of the Oceans 2014 (GEBCO14; Weatherall et al., 2015), which has a horizontal resolution of 30 arc-sec ( ∼ 900 m), referenced to the local MSL. GEBCO data, however, gives a maximum depth of ∼ 300 m in the center of Palu Bay, which is largely in error as this depth is greater than 800 m in more accurate data sets ( Figure 6); hence using GEBCO data may have caused errors in this earlier modeling study. The Indonesian Geospatial Information Agency (BIG) provides a national bathymetric dataset, referred to as BATNAS, which has a horizontal resolution of six arc-sec ( ∼ 180 m) and is referenced to MSL. When comparing BATNAS to geo-referenced satellite images from Google Earth TM , however, the coastline was not accurately located ( Figure 7). BIG also provides bathymetric contours measured in Palu Bay during 2014, 2015, and 2017 surveys (referred to as BIG14, BIG15, and BIG17), all referenced to the lowest astronomical tide. These data sets are shown and discussed in detail by Liu et al. (2020), who point out that BIG17 is the most detailed, but only covers a small portion of the northern section of the bay. They also note that BIG15 is lower resolution and has a few anomalies compared to BIG14, concluding that BIG14 should be used to cover the areas of the bay not covered by BIG17. We proceeded similarly in this work. Regarding topographic data, BIG's national topographic digital elevation model, referred to as DEMNAS, has a horizontal resolution of 0.27 arc-sec ( ∼ 8.3 m) and a vertical datum EGM2008. In contrast to the BATNAS data set, the DEMNAS topography data set is fully consistent with the BIG14 bathymetric contours and satellite images, with both agreeing well at the coast (Figure 7). Finally, as discussed earlier, three post-tsunami bathymetric surveys were reported by Takagi  Thus, in our computational grids, we interpolated the deepwater bathymetry of Frederik et al. (2019) with the shallow water bathymetry from the BIG14 contours, and the topography from DEMNAS, after referencing them all to the same MSL +1 m vertical datum. To avoid numerical instabilities caused by slight discontinuities in bathymetry from combining different datasets, we applied a 2D Gaussian smoothing filter with a standard deviation of 1. The resulting bathymetry and topography are shown in Figure 6.

Numerical Models and Tsunami Modeling Methodology
Landslide tsunamis The 3D non-hydrostatic model NHWAVE is used to simulate initial wave generation and propagation for deforming submarine/subaerial slides, represented by a bottom layer of granular material (debris flow) (Ma et al., 2012;Ma et al., 2015), which was supported by various observations in Palu (e.g. Liu et al., 2020). Euler equations are solved in the water, in a horizontal Cartesian grid (x, y) with boundary fitted σ-layers in the vertical direction and, in the slide layer, conservation equations are depth-integrated. These equations are coupled along the slide-water interface through kinematic and dynamic conditions. One σ-layer achieves the same level of dispersive properties as a 2D Boussinesq model, but more layers allow accurately modeling wave dispersion in larger depth to wavelength ratios. Here we use five layers, which was shown to be adequate for simulating tsunamis from coastal landslides (e.g., Grilli et al., 2015;Schambach et al., 2019). As Palu Bay has steep shores, we use the latest implementation of NHWAVE, in which effects of vertical acceleration are included in the slide layer (Zhang et al., 2021a;Zhang et al., 2021b); these were found important on steep slopes . In the absence of site-specific information, we used the same granular density and internal friction values for the slide material as in Nakata et al. (2020), and a basal friction value at the lower end of their tested range (2-6 deg); this was also the value Grilli et al. (2019) used to model the 2018 Anak Krakatau flank collapse. Hence we have,ρ S 2, 050 kg/m for granular material density, with a 60% solid volume fraction, ρ w 1, 025 kg/m for water density, internal friction angle ϕ i 30 deg, and basal friction angle ϕ b 2 deg. We did not perform a sensitivity study to basal friction, as we did not expect large effects of a small change in friction due to the short distances of slide motion and the rapidly increasing water depth across Palu Bay.
NHWAVE was extensively validated for a variety of tsunami benchmarks (Zhang et al., 2017), including laboratory experiments for slides made of glass beads performed by some of the authors . The model was also used to simulate historical case studies, for which tsunami coastal impact had been measured (Tappin et al., 2014;Grilli et al., 2019;Schambach et al., 2020). In the latter cases, the initial unfailed landslide geometry is first recreated by moving the failed landslide material upslope. The model then simulates both the down-slope motion of the failing slide, coupled to that of the overlying water. For all benchmarks or actual events, NHWAVE was found to perform well and to adequately reproduce the reference data, provided the discretization was sufficient.
In the present simulations, except for slide LS-F* for which we use the actual mapped slide geometry, as in earlier work (e.g., Enet and Grilli, 2007;Grilli et al., 2015;Schambach et al., 2019), the initial slide geometry is modeled at its unfailed location as a sediment mound of quasi-Gaussian crosssections, with maximum thickness T, and an elliptical footprint of down-slope length b and cross-slope width w; with these definitions, the slide volume is calculated as,V S 0.3508 bwT (see Appendix in Schambach et al. (2019), for details). Table 2 Figure 1B). In simulations, the initial geometry of each landslide is carved out of the pre-failed bathymetry BIG14, gridded in NHWAVE's 30 m grid BG, at each slide estimated initial location (Table 2; Figure 6). Since the post-failure coastal bathymetry did not show clear slide deposits, no material was removed from the downslope bathymetry prior to simulations. Based on the shear wave travel time from the earthquake epicenter to Palu Bay discussed above (on the order of 10-15 s), all the landslides should have been affected by ground shaking within seconds of each other; hence, they are all assumed to fail at the same time in the model. However, there was a delay between the first instance of ground shaking due to seismic waves and the actual landslide initiation of motion, likely because complete failure required a sufficient built-up of excess pore pressure (and perhaps some liquefaction) in the submerged toe of the slide material (e.g., Tappin et al., 2008). This delay was estimated to t S x75 s based on the aircraft pilot video and using NHWAVE simulations to compute how much time was required to achieve the observed wave generation. Figure 2 shows the state of wave generation at t 108 s and, modeling the largest slides (in particular LS-F*), we find that it takes 30-35 s of wave generation to qualitatively achieve the same stage as observed in the video, hence on average t S x108 −33 75 s.
On this basis, the generation of landslide tsunamis and their initial propagation up to t t f 150 s were simulated in grid BG with NHWAVE, simultaneously for all the considered slides ( Table 2; an animation of this simulation video4.mp4 is provided in supplementary material). Results show that, at time t f , slides are no longer tsunamigenic and maximum landslide tsunami runups have occurred onshore of each slide location. For t > t f , simulations are continued in FUNWAVE for landslide tsunamis alone or in combination with coseismic tsunamis, based on NHWAVE results for surface elevation and horizontal velocity (interpolated at 0.531 times the local depth for consistency with FUNWAVE). This is detailed next.
Coseismic or dual landslide/coseismic tsunamis Three different coseismic tsunami sources are simulated in grid BG for the 2018 Palu event. Two are modeled for t ≥ 0 (Jamelot et al., 2019;Socquet et al., 2019) with the 2D fully nonlinear and dispersive Boussinesq model FUNWAVE (Wei et al., 1995;Shi et al., 2012), initialized with a static surface elevation equal to the maximum seafloor deformation, and one (Ulrich et al., 2019) using NHWAVE for t ≤ 60 s based on directly specifying the bottom deformation in time and space, and then with FUNWAVE for t > 60 s (see details of coseismic sources later). For the dual earthquake/landslide sources, NHWAVE results are linearly superimposed at t t f with those of FUNWAVE for the simulation of each coseismic source (i.e., surface elevation and horizontal velocity). Simulations of the combined tsunamis are then continued in FUNWAVE for t > t f . FUNWAVE has been extensively validated in various wave propagation and coastal inundation studies (e.g., Shi et al., 2012), including for tsunami inundation/runup and coastal velocity benchmarks (Horrillo et al., 2015;Lynett et al., 2017), and applied to tsunami case studies, both historical and hypothetical (e.g. Tappin et al., 2014;Grilli et al., 2015;Grilli et al., 2019;Schambach et al., 2019;Schambach et al., 2020).
For all cases simulated here, landslide or coseismic tsunamis alone, or dual sources, FUNWAVE simulations are performed by one-way coupling in the 2-level nested Cartesian grids ( Figure 6) BG (30 m resolution) and EG/SG (7.5 m resolution; see the earlier studies for details). To prevent reflection at open boundaries, 1800/4,200 m wide sponge layers are specified along the western/ northern boundaries of grid BG. Inundation and runup are modeled along coastal boundaries by way of a moving shoreline algorithm, with dissipation by wave breaking and bottom friction being simulated in FUNWAVE ; here, a bottom friction coefficient C d 0.0025 is used, which corresponds to coarse sand (also used in NHWAVE).
As discussed in the introduction, all published studies of 2018 Palu to date used non-dispersive tsunami propagation models (Table 1). However, earlier work has shown the importance of frequency dispersion when modeling landslide tsunami generation and propagation, particularly when the initial slide footprint, and hence initial wavelength, are small compared to depth (e.g., Ma et al., 2012;Schambach et al., 2019). Here, we use dispersive models (NHWAVE and FUNWAVE) and, to assess the importance of dispersive effects, we run some simulations by turning off dispersion terms in the models (results are detailed later). In Palu Bay, while landslide tsunami waves generated in very shallow water may not initially be significantly dispersive, dispersion would become larger once waves propagated into the much deeper water of the center of the bay. Dispersion affects both phase speed and wave-wave interactions during propagation and, ultimately, tsunami coastal impact. Additionally, close to shore, dispersion may create undular bores (a.k.a. dispersive shock waves) near the crest and trough of longer shoaling tsunami waves, enhancing coastal impact (e.g., Madsen et al., 2008). This was demonstrated by Schambach et al. (2019) who also showed that higher resolution nearshore grids must be used to capture undular bores in FUNWAVE, which is one of the reasons here for using the 7.5 m grids EG and SG, even though the bathymetric data was not available at that level of detail; but, wave physics may call for it.

Earthquake Source Models
As there is no consensus on the 2018 Palu earthquake parameters, which fault(s) was(were) responsible, and how the rupture proceeded, to assess the source-related epistemic uncertainty in tsunami simulations, we modeled three representative earthquake sources, whose main characteristics are summarized in Table 1.
The first two sources use multiple subfaults whose parameters (depth, dimension, angles) were optimized, using Okada (1985)'s method, to best match observed onland displacements inferred from pre-and post-earthquake satellite observations, while satisfying the M7.5 centroid moment tensor. Jamelot et al. (2019) thus define nine subfaults with specified length and strike, and use displacements derived from Sentinel-2 satellite images to optimize other parameters. Socquet et al. (2019) use 294 subfaults (42 in the strike direction and seven in the dip direction), and displacements inferred from Sentinel-2 and Landsat-8 satellite images, as well as SAR interferograms from ALOS-2 satellite data to optimize fault parameters. For these two sources, we computed the maximum seafloor deformation with Okada (1985)'s model and used it as initial surface elevation in FUNWAVE (with no initial velocity). As some studies indicated that effects of steep shores might have been important (Carvajal et al., 2019;Heidarzadeh et al., 2019;Jamelot et al., 2019;Ulrich et al., 2019), we computed the additional vertical displacements due to horizontal displacements using Tanioka and Satake (1996)'s method. Assuming a supershear rupture, we specified the initial time for these surface elevations in FUNWAVE to t 15 s into the event. Figures 8A,B show the initial elevations computed for these sources, which clearly are aligned along different fault traces (Figure 1), but both predict a positive initial elevation (seafloor uplift) on the east side of the bay and a negative initial elevation (seabed subsidence) on the west side. Jamelot et al. (2019)'s source was designed by the authors to cause larger elevations in the area of the Pantoloan tide gauge and Wani to the north of the bay, and in the area of Grand Mall and the KN Hotel to the south in Palu City, where large tsunami impact was observed. This is clearly reflected in Figure 8A, by the larger initial elevations for this source in those areas.
The third source by Ulrich et al. (2019) was developed from physics-based modeling of the earthquake failure in space and time using Seisol, which solves elastodynamic wave equations for spontaneous dynamic ruptures and seismic wave propagation (Dumbser and Käser, 2006;Pelties et al., 2014;Uphoff et al., 2017). Model inputs included geometry and frictional strength of the fault, tectonic stress state, and regional lithological structure, which were determined from datasets specific to the Palu region. Based on the authors' results for the horizontal and vertical ground motions (provided every 0.5 s for 50 s over a dense grid), we created time-series of seabed motion, which were used as bottom boundary conditions over grid BG in NHWAVE. As with the other sources, contributions to vertical displacement due to the horizontal movement of steep slopes were included. Simulations of tsunami generation/propagation were done in 3D with NHWAVE up to t 60 s and then in 2D with FUNWAVE for t > 60 s. To compare results with the other two sources, the surface elevation computed with NHWAVE is plotted at t 15 s in Figure 8C, where, we see that, while the deformation is aligned along a fault trace similar to that of Socquet et al. (2019), likely due to their very physics-based modeling, very different from that of the other two sources, the polarity of deformation is opposite, i.e., there are large negative elevations (subsidence) on the east and large positive elevations to the SW and NE, of the bay.

RESULTS
Simulations were performed with NHWAVE and FUNWAVE following the methodology detailed above, by considering first each of the three coseismic sources (Figure 8), then the seven parameterized landslide sources ( Table 2), and finally for dual earthquake/landslide sources combining each coseismic source with all the landslide sources. Simulations for earthquake or landslide sources alone were performed with or without dispersive effects in the models. All simulations were performed up to t 1, 200 s, which was determined to be long enough for maximum runup to be achieved along the Palu Bay shores. Similar results were computed for each type of simulation: 1) the envelope of maximum surface elevations over the study area and runups along the Palu Bay coastline, to be compared with measurements from field surveys (Pribadi et al., 2018;Mikami et al., 2019;Omira et al., 2019;Putra et al., 2019;Widiyanto et al., 2019), in Figures 3,9,11, for the coseismic, landslide, and dual sources, respectively (runups computed in both the coarser BG and finer SG/EG grids are provided, whereas envelopes are computed using the finer resolution results, wherever available); and 2) time series of surface elevations computed at locations where various observations or measurements were made or inferred, i.e., Wani, Pantoloan, Dupa, Talise, KN Hotel, and Grand Mall ( Figure 1B), in Figures 4,10,12, respectively.
First, regarding the effects of grid resolution, Figures 3,9,11 show that for all types of sources the largest runups simulated in the finer grids (EG/SG) are larger than in the coarser grid (BG), which justifies using nested grids in FUNWAVE. Then, regarding dispersion, for coseismic tsunamis, time series of surface elevation in Figure 4 show that dispersion causes only small absolute changes in wave elevation (mostly near the crests), at most times and locations, compared to non-dispersive simulations. This is expected for the longer coseismic tsunami waves; nevertheless, relative differences between the two simulations can be 25-35% at some times/locations. In contrast, for landslide tsunamis, Figure 10 shows that dispersion causes much larger absolute or relative differences in surface elevations, and larger phase shifts. This is also expected, based on earlier work (e.g., Glimsdal et al., 2013;Schambach et al., 2019), for shorter landslide tsunami waves, particularly considering the large depth of the bay. These results justify using dispersive wave models in this work.
For coseismic sources alone, Figure 3 shows a similar trend for runups predicted around Palu Bay, but with large absolute differences; in particular runups are in general lower for the Socquet et al. (2019) source. All three sources, however, significantly underpredict runups observed in the southern part of the bay (south of ∼ − 0.75 + N), with a maximum of 5 m whereas observed runups reached up to 10.5 m. In contrast, runups are relatively well predicted in the northern part of the Bay by Jamelot et al. (2019)'s andUlrich et al. (2019)'s coseismic sources, with a slight advantage for the latter. This agrees with conclusions of earlier studies that coseismic sources alone cannot explain the tsunami coastal impact, particularly in the south (e.g., Nakata et al., 2020;Sepúlveda et al., 2020). Consistent with these results, Figure 4 shows that measured or inferred time series of surface elevation are not well reproduced, particularly in the southern part of the bay. Exceptions are Wani, the northern location ( Figure 1B Figure 3 differ from those published by the sources' authors or others who used these. Reasons for these differences are multiple: 1) we use higher resolution model grids based on higherresolution bathymetric data than in the earlier studies, hence both wave physics and runup are more accurately simulated; 2) unlike earlier studies, we use dispersive wave models in which wave-wave interactions differ, even for coseismic sources; 3) for Ulrich et al. (2019)'s source, we generate the tsunami dynamically for 60 s in a 3D model, whereas they used a 2D NSWE model; and 4) finally FUNWAVE has a particularly accurate moving shoreline algorithm to capture runup on steep slopes, which may not be the case for all models.
For the landslide sources alone, Figure 9 shows that observed runups are well predicted in the SW part of Palu Bay, particularly in the area of the largest landslide sources (LS-E, LS-F*). However, a few of the largest observed runups are still underpredicted in the area of Dupa on the SE of the Bay (around −0.85 + N). In the northern part of the bay, on the western side, observed runups are nearly as well predicted as for coseismic sources, but because of the timing of the event,  maximum runups caused by coseismic or landslide sources would have likely occurred at different times here (see animations of model results). Hence, it is difficult to identify their primary source, which may explain the mitigated conclusions or even confusion in some earlier studies. On the NE side of the bay, around Wani and Pantoloan, the landslide tsunami impact is predicted to be quite large and explains the large runups observed better than for coseismic sources. This is confirmed in time series of surface elevation (Figure 10), where there is a much better agreement in Wani and Pantaloan of model results with the reconstituted time series than for the coseismic sources. In the SE of the bay, however, consistent with the underpredicted runups, the landslide tsunami simulations do not explain well the time series reconstructed in Dupa, Talise, and KN Hotel. Finally, at Grand Mall in the south, while the shorter period landslide tsunami waves agree better in timing with those of the reconstructed time series, their amplitude is still underpredicted, despite using a very fine model grid that could have enhanced wave shoaling.
Results of the dual earthquake/landslide source simulations ( Figure 11) are consistent with the above observations. In all cases, but particularly for the combination of Ulrich et al. (2019)'s with the landslide sources, the observed runups on the entire west side of the bay are well simulated at most locations, and this is also the case on the east side of the bay, except for the area around Dupa where another local source of waves is required, perhaps from another landslide not yet identified in this region, as pointed out in some other studies (e.g., Liu et al., 2020). In the time series results of Figure 12, we see a good agreement between the dual source simulations with the reconstructed time series at Wani and Pantoloan, particularly when using Ulrich et al.  Table 2) sources, compared with field measurements (see Figure 3 for definitions of data, three coseismic sources, line colors and types). Landslide tsunami generation is first simulated with NHWAVE in grid BG (figure footprint) up to t f 150s. At this time, NHWAVE results are linearly combined (surface elevation and velocity at 0.531 times the local depth) with those computed with FUNWAVE in grid BG for each of the three coseismic sources. Simulations are then continued with FUNWAVE for t > t f , in grid BG and then nested EG/SG grids (white footprints in figures (B,D,E)).

DISCUSSION
Our work here shows that for the 2018 Palu event, a combination of earthquake and coastal landslides generated the tsunami. We also show that mapped (rather than hypothetical) landslides were critical in achieving this result. Video evidence was also instrumental in differentiating between the two possible mechanisms, especially at Pantoloan, where the tide gauge data used to validate previously published numerical models, was found to be partly misleading because it filtered out the high frequency landslide tsunami waves. To confirm this, we applied to the model results of Figure 12B a 30 s moving average and 1 min downsampling similar to that of the tide gauge Frontiers in Earth Science | www.frontiersin.org January 2021 | Volume 8 | Article 598839 17 (Sepúlveda et al., 2020). Supplementary Figure S2 (in supplementary Data Sheet 1) shows that this eliminates shorter waves from the time series, such as caused by the landslides or seen in the video recording near Pantoloan dock. Additionally, the filtered results based on Ulrich et al.
(2019)'s dual source agree well with the first few waves in the tide gauge record, although amplitudes are smaller.
Our simulations of published earthquake sources (Jamelot et al., 2019;Socquet et al., 2019;Ulrich et al., 2019) show the epistemic uncertainty associated with modeling the coseismic tsunami. While the initial surface elevation from each coseismic source is quite different (Figure 8), the generated tsunamis all reproduce the runups observed in the northern section of the bay, but underpredict the larger runups in the south. Without a comparison of pre-and post-earthquake leveling data, it is not clear which, if any, of the earthquake models is most appropriate. The recent work by Natawidjaja et al. (2020), published too late to include for consideration here, reinterpreted Frederik et al. (2019)'s multibeam bathymetry and identified the major, meandering, submarine channel in the center of Palu Bay, as the seabed expression of the 2018 movement of the strike-slip fault. Based on this study, the fault could be considered to be more effective in tsunami generation than previously proposed. Several aspects of their model, however, lead us to conclude that further justification is required before it can be accepted as a viable alternative to those already published: 1) it is so very different to previously published interpretations based on the same datasets (Frederik et al., 2019;Liu et al., 2020); 2) interpretations of several meters of vertical seabed movement in the context of the resolution of Frederik et al. (2019)'s bathymetry ( ± 5 m) is questionable, as is the identification of recent (2018) seabed movement; and 3) the suggestion that the earthquake triggered "massive" SMFs to account for the tsunamis is contradicted by the available bathymetric evidence (Frederik et al., 2019;Liu et al., 2020).
In the southwest, the large runups observed just onshore of confirmed coastal landslides stress the importance of simulating landslide tsunamis. In this context, using a dispersive numerical tsunami model was particularly important for accurately propagating the shorter wavelength landslide tsunami waves. Here, Liu et al. (2020)'s shallow water bathymetric survey was critical in parameterizing the numerous coastal landslides. The detailed surveying of slide LS-F* by Takagi et al. (2019), where large wave generation was observed (see Figure 2), provided additional information. These authors were the first to identify these landslide tsunami mechanisms upon which we built our more complex and comprehensive model. The videos were instrumental in allowing us to identify a landslide trigger delay of 75 s, with the pilot video and time series of surface elevation at Wani and Pantoloan ( Figures 10A,B) providing key evidence. The mapped landslildes we used in our modeling are found to be capable of generating runups on the same order as those observed onshore of their locations and their reconstructed time series impact, with good agreement at most locations. One exception is in the southeast of the bay, where runups are still underpredicted in the Dupa area (simulated 2-4 m, vs observed 8-10.5 m). At the Grand Mall, wave arrival matches that observed, however, the amplitudes are not as large.
To explain the large runups observed in the SE of the bay, Nakata et al. (2020) modeled a large 700 ×10 6 m hypothetical SMF off of Talise and obtained a good agreement with observations near this location. There is no indication on the seafloor for such a large recent failure, although Liu et al. (2020) suggest that there are several large SMFs south of this location, but these do not appear to be recent. Our simulations suggest that an additional SMF off of Talise could explain both the large runups observed between Dupa and Talise, and improve the agreement of simulations with the time series inferred at Talise and the KN Hotel. However, as our stated goal was to only model proven landslide sources, we did not consider such a hypothetical SMF in our earlier dual sources.
Nevertheless, to test this hypothesis, we modeled a SMF at Nakata et al. (2020)'s location (i.e., 119.8675+ Lon. E, −0.8540+ Lat. N), where Liu et al. (2020) identify recent seabed movement, albeit with a much smaller volume. This SMF's elliptical footprint is marked in Figure 13A, with dimensions b 500 m by w 1,000 m; given a maximum thickness T 150 m, the SMF volume is,V s 26.3 ×10 6 m. As with the other landslides, the tsunami was first simulated with NHWAVE without a coseismic source, and then propagated with FUNWAVE in grids BG and SG. Figures 13A,C show that the landslide tsunami focuses on two areas onshore of the SMF: 1) just south of the KN Hotel (−0.866+ N) causing a 3 m runup, consistent with measurements; and 2) north of Talise Table 2. Some wave interferences slightly reduce the runup south of the KN Hotel, but north of Talise, the combined runup is still nearly 6 m, whereas without the SMF it was only 3 m ( Figure 11C). Finally, Figures 13E,F show time series of surface elevations computed at Talise, KN Hotel, and Grand Mall, respectively. Compared to earlier results in Figure 12, the new dual source simulations improve the overall agreement with reconstructed time series. At the KN Hotel, in particular, only the inclusion of the SE SMF can explain the leading elevation wave observed at t 125 s (underestimated but arriving at the correct time).
The 2018 Palu tsunami was unusual, and complicated, with a strike-slip earthquake mechanism which triggered coastal landslides. Previous publications show how difficult it has been to identify the tsunami generation mechanism(s). Our work here, however, demonstrates that, for most of Palu Bay, the earthquake and the mapped coastal landslides were equal contributors to the large runups measured around the bay, except in the southeast where an additional (although partly hypothetical) SMF is required. We show the importance of modeling dual earthquake/landslide sources, and of considering all available information to identify how the tsunami waves were generated including, for the first time, time series of tsunami impact reconstructed from video evidence, in addition to the (normally used) runup evidence from field surveys, tide gauge data, and survivor accounts.
A proper understanding and modeling of such destructive dual source tsunami events can help mitigate tsunami coastal hazard resulting from future similar events, here or in other tsunamiprone areas.

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: https://drive.google. com/drive/folders/1kPUnYcenRFa0KLhzhyZgJG9eya5CfC8D? usp sharing.

AUTHOR CONTRIBUTIONS
LS performed all of the tsunami simulations, including postprocessing, and worked on the manuscript. SG supervised all aspects of the work and worked on the manuscript. DT provided insight into the marine geology and worked on the manuscript. resources, as part of the Extreme Science and Engineering Discovery Environment (XSEDE) (project BCS-170006), which is supported by the National Science Foundation (NSF) grant number ACI-1548562.

ACKNOWLEDGMENTS
Special thanks to Anne Socquet for providing their earthquake slip distribution model, and to Thomas Ulrich for providing the time series of horizontal and vertical ground motion from their earthquake model. DT publishes with the permission of the CEO of the British Geological Survey United Kingdom Research and Innovation. Boma Kresning is acknowledged for helping to navigate and translate Indonesian data sources. The Pantoloan tide gauge records and related informations were obtained from the Agency for Geo-spatial Information, Indonesia (BIG) (http://tides.big.go.id).