Insights on the Source of the 28 September 2018 Sulawesi Tsunami, Indonesia Based on Spectral Analyses and Numerical Simulations

The 28 September 2018 Sulawesi tsunami has been a puzzle because extreme deadly tsunami waves were generated following an Mw 7.5 strike-slip earthquake, while such earthquakes are not usually considered to produce large tsunamis. Here, we obtained, processed and analyzed two sea level records of the tsunami in the near-field (Pantoloan located inside the Palu Bay) and far-field (Mamuju located outside the Palu Bay) and conducted numerical simulations to shed light on the tsunami source. The two tide gauges recorded maximum tsunami trough-to-crest heights of 380 and 24 cm, respectively, with respective dominating wave periods of 3.6−4.4 and 10 min, and respective high-energy wave duration of 5.5 and >14 h. The two observed waveforms were significantly different with wave amplitude and period ratios of ~16 and ~3, respectively. We infer tsunamigenic source dimensions of 3.4–4.1 km and 32.5 km, for inside and outside of the Palu Bay, respectively. Our numerical simulations fairly well reproduced both tsunami observations in Pantoloan and Mamuju; except for the arrival time in Mamuju. However, it was incapable of reproducing the maximum reported coastal amplitudes of 6–11 m. It is possible that these two sources are different parts of the same tectonic source. A bay oscillation mode of ~85 min was revealed for the Palu Bay through numerical modeling. Actual sea surface disturbances and landslide-generated waves were captured by two video recordings from inside the Palu Bay shortly after the earthquake. It is possible that a large submarine landslide contributed to and intensified the Sulawesi tsunami. We identify the southern part of the Palu Bay, around the latitude of -0.82oS, as the most likely location of a potential landslide based on our backward tsunami ray tracing analysis. However, marine geological data from the Palu Bay are required to confirm such hypothesis.


Introduction
The Sulawesi Island of Indonesia was the site of a large earthquake (Mw 7.5) on 28 September 2018 which was followed by a deadly tsunami. The United States Geological Survey (USGS) located the epicenter at 0.178 o S and 119.840 o E occurring at 10:02:43 UTC at the depth of 10.0 km (Fig. 1). The USGS focal mechanism solution identified a dominant strike-slip mechanism with strike, dip and rake angles of the most likely fault plane as 350 o , 67 o and -17 o , respectively (red mechanism in Fig. 1). The fault plane from the gCMT (global CMT: http:// www.globalcmt.org) were: 348 o , 40 o and -9 o , respectively. Based on media reports and results of preliminary field surveys (Figure 2), the Sulawesi tsunami reached a maximum height of approximately 6-11 m along the coast of Palu within the Palu Bay ( Fig. 1)  ). As of 10 October 2018, the reported number of deaths has been 1763, making the Sulawesi event as the deadliest event in Indonesia since the December 2004 Sumatra-Andaman earthquake and tsunami. It is likely that the number of deaths will increase in the following weeks. Indonesia has experienced several deadly tsunamis since the December 2004 tragedy ([ 130,000 deaths;Synolakis and Kong, 2006) including the tsunamis of October 2010 Mentawai (408 deaths) (Satake et al. 2013), the 2006 West Java (668 deaths) (Fujii and Satake, 2006), and the March 2005 Nias-Simeulu ([ 915 deaths, mostly from earthquake) (Borrero et al. 2011).
From a tectonic point of view, the epicentral area around the Sulawesi Island is located within one of the most complicated tectonic zones of the world where four major tectonic plates meet: Pacific plate from the east, the Philippine Sea plate from north, the Australian plate from south and the Sunda plate from the east (Fig. 1). The Palu-Koro strike-slip fault system (PKF) runs through the epicentral area (Fig. 1a). The Sulawesi region experienced at least three other lethal tsunamis in past hundred years (  The 28 September 2018 Sulawesi event has been considered as a tsunami mystery because strike-slip earthquakes are usually unable to produce large tsunamis (Synolakis, 2003; while the 2018 Sulawesi tsunami was 6-11 m in height along the coast. Strike-slip ruptures do not usually produce large vertical uplift/subsidence on the seafloor because of the dominant horizontal movements of the fault; thus, they are usually unable to produce destructive tsunamis. For example, the Mw 7.8 Wharton Basin earthquake in SW Sumatra, with a dominant strike-slip mechanism, generated tsunami amplitudes of less than 10 cm at tide gauges ). In addition, the size of the earthquake was Mw 7.5 which is relatively small in terms of tsunami generation. Therefore, it has been speculated that the large waves of the 2018 Sulawesi tsunami were the results of secondary sources triggered by the main shock. Potential secondary sources are splay faulting, submarine landslides, bay oscillations and funneling effects of the Palu Bay ( Fig. 1).
In this research, we present the observed waveforms of the tsunami and conduct analysis to explore the spectral contents of the waveforms. Numerical modeling of the tsunami is performed using the tectonic source of the tsunami in order to examine whether the tsunami could be generated by the purely-tectonic source or not. Results from spectral analysis are combined with numerical simulations to infer potential sources of the tsunami.

The Earthquake Fault and Local Bathymetry
The Palu-Koro strike-slip fault system (PKF) was most likely responsible for the Mw 7.5 Sulawesi earthquake (Fig. 1, 3). According to Bellier et al. (2001), the left-lateral and north-south striking PKF   (Fig. 3a). We note GEBCO bathymetry data may not be accurate enough for coastal and shallow areas such as the Palu Bay. The bathymetry data from BIG (Badan Informasi Geospasial, Indonesia) indicates a maximum water depth of * 700 m (Fig. 3b) as a deep narrow channel in the middle of the bay.

Data and Methods
We analyzed three tide gauge records in Mamuju (Indonesia), Pantoloan (Indonesia) and Lahat Datu (Malaysia). All records have sampling intervals of 1 min. The Lahat Datu data was provided by the Intergovernmental Oceanographic Commission of UNESCO/UN at http://www.ioc-sealevelmonitoring. org/ while the other two records were obtained and processed by the Agency for Geo-spatial Information, Indonesia (BIG) (http://tides.big.go.id). Sea level data were de-tided by predicting the tidal signals using a least squares method of harmonic analysis and then removing the tidal signals from the original tide gauge records. For spectral analysis, we followed the Welch algorithm (Welch, 1967) with half-window overlaps (e.g., Rabinovich et al. 2011;Heidarzadeh and Satake, 2013). Duration of high-energy tsunami waves were calculated by the Averaged-Root-Mean- Square (ARMS) diagrams for the tsunami waveforms (Heidarzadeh and Satake, 2014) using a moving time window with length of 20 min (i.e., 20 data points). The duration of high-energy tsunami wave is the time that the ARMS level of the waveform is above that for the background wave (i.e., before the arrival of the tsunami).
For tsunami simulations, we applied the COM-COT package (Cornell Multi-grid Coupled Tsunami Model) (Liu et al. 1998;Wang and Liu, 2006). COMCOT solves linear and nonlinear Shallow Water Equations on both Cartesian and Spherical coordinates. COMCOT adapts a numerical scheme to mimic physical dispersion following the method by Imamura et al. (1988). A two-level nested bathymetric grid was created with resolutions of 30 and 5 arc-sec from the GEBCO-2014 (The General Bathymetric Chart of the Oceans) data (Weatherall et al. 2015) (Fig. 1a). The higher-resolution grid (i.e., level 2) was interpolated from the original GEBCO 30 arcsec bathymetry and topography data. We applied linear simulations for the larger grid and nonlinear simulations for the smaller ones. Runup calculations were not included in our numerical simulations, as it requires high-resolution nearshore bathymetry and topography. Simulations were conducted for a total time of 11 h with time step of 0.5 s. The second version of the USGS source model has been used in this study for tsunami generation which has dimensions of 200 km (length) 9 30 km (width) with maximum slip of 9 m occurring within the Palu Bay (https://earthquake.usgs.gov/earthquakes/eventpage/ us1000h3p4/finite-fault). The initial tsunami wave field was calculated using the co-seismic seafloor displacements obtained from the Okada's (1985) analytical method. Both vertical co-seismic displacement (D z ) and the vertical water displacement due to steep slopes (D h ) were considered following the approach proposed by Tanioka and Satake (1996): in which, D x and D y are the E-W and N-S components of co-seismic seafloor displacement, H is water depth (positive downward), D h is the vertical displacement of water due to the horizontal movement of the slope. The combined vertical displacement of water is the sum of D z and D h .

Results of Tide Gauge Data Analysis and Physical Properties of the Tsunami
The Lahat Datu tide gauge did not record any clear tsunami signal (Fig. 1b) which can be attributed to the relatively large noise level (± 8 cm) and small size of the tsunami. The maximum trough-to-crest wave height recorded at two stations of Pantolaon (inside the Palu Bay) and Mamuju (outside the Palu Bay) are 380 and 24 cm, respectively (Table 1, Figure 1c). The Pantoloan record is 16 times larger than that of the Mamuju. The respective arrival times are 5 and 19 min (Table 1). Our numerical simulations, however, indicate an arrival time of *45 min for the Mamuju record. Therefore, either the timing of the tide gauge was problematic or maybe a different mechanism was responsible for the early tsunami observed in Mamuju. Assuming that the latter holds true, the relatively short tsunami arrival time in Mamuju may indicate that other sources nearby were responsible for the tsunami in Mamuju.
In addition to such a significant difference in tsunami heights between Pantoloan and Mamuju, the two records are significantly different in terms of wave periods: their respective dominant periods are 3-4 and 10-12 min from inspection of the waveforms (Fig. 1c, Table 1). The dominant period of the Pantoloan waves are *3 times longer than those of Mamuju. The sampling interval of 1 min appears to be insufficient to properly capture the short-period tsunami oscillations in Pantoloan; thus, the actual extrema of the waves are possibly higher than the values shown in Table 1. The ARMS diagrams (Fig. 4) show that the high-energy waves at the Pantoloan last for 5.5 h while they last for at least 14 h in Mamuju.
Tide gauge data analysis indicates that, in terms of physical properties of the waves, the Pantoloan and Mamuju tsunami waves are significantly different: the Pantoloan waves are significantly larger in height and are dominated by shorter-period waves and shorter duration in comparison to the waves at the Mamuju station. It appears that the two waveforms are possibly generated by two completely different sources or mechanisms. Usually different tide gauge records of the same tsunami contain the same or   Tanioka and Satake (1996) Vol. 176, (2019) Insights on the source of the 28 September 2018 Sulawesi tsunami, Indonesia 31 similar wave periods for at least the first few cycles (e.g., Heidarzadeh and Satake, 2013).

Results of Spectral Analysis
Spectral analyses (Fig. 5) reveal dominant tsunami period bands. For the Pantolaon station, the tsunami period band is 2.8-6.7 min while it is 4.0-18.3 min for Mamuju. The respective dominant tsunami periods are 3.6-4.4 and 10.0 min (Fig. 5). Tsunami source dimensions can be estimated from the tsunami dominant periods using the following equation (Heidarzadeh and Satake 2015b): in which, L is source length, H is average water depth, g is gravitational acceleration and T is tsunami dominant periods. We note that Equation (2) is applicable for wave periods (T) extracted from first few cycles of the tsunami, which is usually free from bay oscillations and other bathymetric effects. By assuming average water depth of *100 m for inside the Palu Bay (Fig. 2) and *1200 m for outside the bay (Fig. 1), Eq.
(2) results in tsunami source dimensions of approximately 3-4 and 32-38 km for inside and outside the Palu Bay, respectively. We note that the average water depth of the Palu Bay could be slightly different considering the BIG bathymetry data in Fig. 3b. It is possible that these two sources are different parts of the same tectonic source because the isolated segment of the tectonic source located inside the Palu Bay may act independently towards tsunami generation within the bay.

Results of Numerical Simulations
Results of numerical simulations are shown in Figs. 6-11. Two different initial conditions for the 2018 Sulawesi tsunami were considered: the coseismic vertical crustal displacement (Fig. 6c, D z ), and the combination of the vertical component (D z ) and the additional vertical water displacement due to the horizontal movement of the seafloor slope (D h ) (Fig. 6d, the component D z þ D h ). It can be seen that these two initial conditions are similar to each other with small differences (Fig. 6). Our calculations yield up to *3 and *1 m of co-seismic horizontal N-S (Fig. 6a) and E-W (Fig. 6b) crustal displacement. The maximum co-seismic vertical displacement is *1 m (Fig. 6c, d).
Snapshots of tsunami simulations at different times for two large and small grids are shown in Figs. 7, 8 indicating a similar pattern of tsunami propagation from our two initial sources (i.e., D z and D z þ D h Þ. Several nodes (i.e., points with zero amplitude) appear in the tsunami snapshots within the Palu Bay (Figs 7, 8, right panels) indicating patterns of bay oscillations. Based on the snapshots in grid-1 (i.e., the largest grid), the first negative wave arrives in Mamuju *45 min after the earthquake origin time while the observation waves show an arrival time of *20 min (Fig. 1c). This may indicate that the timing of the Mamuju tide gauge station is not correct or possibly a different local source was responsible for the early tsunami in Mamuju. As the periods and amplitudes of the simulated and observed waves agree well at the Mamuju station, it is most likely that there was a clock error with the tide gauge. The simulated waves at the locations of the two tide gauges agree well with those of observations ( Fig. 9, blue waveforms), at least for the first cycle of the waves. We note that the observation waveform in Mamuju is delayed 45 min in Fig. 9 to match the timing of the simulations. Tsunami simulations are usually aimed at reproducing the first wave cycle because later cycles are mixed with reflected/refracted waves as well as nonlinear bathymetry effects; hence, high-resolution bathymetry/topography data along with dispersive modeling would be necessary to fully reproduce late cycles. In such a context, the simulation results shown in Fig. 9 are successful in reproducing most of the features of the two observed tide gauge records.
Oscillations of tsunami waves at various locations within the Palu Bay are plotted in Fig. 10  our numerical simulations. A long period of *85 min is seen in all of the numerical wave gauges within the Palu Bay which can be attributed to the fundamental periods of the bay. According to Rabinovich (2010), the fundamental periods of a semi-enclosed long basin (i.e., length is significantly larger than the width) is given by the following equation: in which, T is the fundamental mode of the bay, L is the length of the bay, H is average water depth, g is gravitational acceleration and n ¼ 0; 1; 2; . . . which gives various fundamental modes. By assuming a length of 35 km (L) and an average water depth of *100 m (H), the few first fundamental periods are 75, 25, 15 and 10.7 min. The longest fundamental period of 75 min is close to the value of *85 min seen on the simulated waveforms in Fig. 10. The maximum tsunami wave amplitudes during the entire time of tsunami simulations are shown in Fig. 11 for both initial conditions. Figure 11, which does not include runup stage, yields maximum amplitude of *1.5 m for the Sulawesi tsunami. Previous studies have shown that maximum coastal wave amplitudes are rough estimates of wave runup. Obviously, this value is much smaller than the reported runup of 6-11 m by tsunami survey teams . Based on Fig. 11, maximum wave amplitudes are concentrated at the southern tip of the Palu Bay as well as around the Pantoloan tide gauge station.

Insights on the Source Mechanism and Potential Landslide Location
Numerical modeling and spectral analysis conducted in this study help to comment on the type of the tsunami source. Based on the spectral analysis, the tide gauge stations located inside and outside of the Palu Bay generated tsunami waves with significantly-different dominating periods of 4.0 and 10.0 min, respectively. Numerical modeling proved that a purely-tectonic source of the tsunami is able to generate both types of the shorter-and longer-period waves inside and outside of the bay. However, the simulated maximum tsunami amplitudes (1.5 m) are smaller than the surveyed runup heights of 6-11 m inside of the Palu Bay. Therefore, it may be concluded that an additional secondary source has contributed to the tectonic source of the tsunami. We provide evidence in the next section that such a secondary source most likely was a submarine landslide within the Palu Bay.
To shed light on the location of a potential secondary landslide source, a backward tsunami ray tracing is conducted whose result is shown in Fig. 12. For this analysis, a hypothetical tsunami source is placed at the location of the Pantoloan tide gauge station and propagation of the tsunami is modelled up to 5 min, which is the arrival time of the 2018 Sulawesi tsunami at this tide gauge station (Table 1). The edges of tsunami travel distances after 5 min in both directions (north and south) determine potential landslide locations (i.e., the two dashed boxes in the bottom-right panel in Fig. 12). The field survey by Muhari et al. (2018) revealed that the largest runup heights were concentrated at the southern part of the Palu Bay around the latitude of -0.82 o S. Therefore, the bottom dashed rectangle at the latitude of -0.82 o S (Fig. 12) is more likely to be the location of a potential landslide.

Discussion
Our purely-tectonic source model for the 2018 Sulawesi tsunami reproduces fairly well both tide gauge records in Pantoloan and Mamuju; except for the arrival time of the Mamuju record. The agreement Figure 11 a Maximum tsunami amplitudes after 11 h of tsunami simulations based on the vertical co-seismic displacement (D z ). b same as ''a'' but for simulations based on the combination of the vertical co-seismic displacement and the vertical water displacement due to steep bathymetry Vol. 176, (2019) Insights on the source of the 28 September 2018 Sulawesi tsunami, Indonesia 37 between our simulations with observations includes both wave periods and wave amplitudes. In other words, the purely-tectonic source model explains the sharp difference between wave amplitudes (i.e., a ratio of *16) and wave periods (i.e., a ratio of *3) of the Pantoloan and Mamuju tide gauge records. However, the maximum simulated coastal amplitudes from our source model is *1.5 m which are far less than the reported wave heights of 6-11 m. As of 20 October 2018, preliminary tsunami field surveys have been conducted and the data on the runup are yet to be analyzed. It has been reported that the extreme wave height of 11 m was the splash height and perhaps the flow depth on the land was *3 m ). If we assume that these speculations hold true, then our source model could be capable of reproducing runup heights of * 3 m given highresolution bathymetry/topography. Long-wave runup significantly depends on the actual nearshore bathymetry/topography, and particularly on the last coastal slope on which runup occurs (Kanoglu and Synolakis 1998).
In early days following the Sulawesi disaster, the potential involvement of a secondary source such as a submarine landslide (e.g., Synolakis et al. 2002) or a splay fault (e.g., Plafker, 1972) was mentioned. Both phenomena have been responsible for extremely large tsunamis following moderate earthquakes in the past. The 15-m large tsunami waves in July 1998 in Papua New Guinea were generated by a landslide tsunami following an Mw 7.0 earthquake (Tappin et al. 2001;Synolakis et al. 2002;Heidarzadeh and Satake, 2015a). Splay fault have been reported to double tsunami heights compared to cases without a splay component (Heidarzadeh et al. 2009). Both splay faults and landslides are local phenomena with limited source dimensions (e.g. a few kilometers). The dimensions of the 1998 PNG landslide were * 4 km in both length and width. For the 2018 Sulawsi event, a video recording of the sea surface of the Palu Bay shortly after the earthquake, captured onboard an aircraft, shows several local sea surface disturbances (Fig. 13). Another 4-min eyewitness video from onboard a boat also recorded at least three subaerial landslides within the Palu Bay following the earthquake (Fig. 14). Due to the small size of these sea surface disturbances and landslides in the eyewitness videos (Figs. 13, 14), they may not necessarily be of the sizes capable of producing significant tsunamis. It appears that the area is susceptible to subaerial and submarine landslides; therefore, it could be possible that a large submarine landslide along the submerged slopes was responsible for the 6-11 m tsunami waves. Potential contribution of such secondary sources (i.e., co-seismic landslides) could be examined by acquisition and interpretations of marine geological data from the seafloor within the Palu Bay.
Ideas of the existence of a tsunami generation mechanism other than the strike-slip fault might be strengthened by findings from field surveys  which revealed that maximum inundation of the tsunami was limited to the areas within 400 m from the shoreline. Due to their short-period waves, landslide-generated waves are not capable of producing long inundations. From the hydrodynamics point of view, inundation distance strongly depends on the period of the waves: the longer the period of the waves, the longer the inundation distance. The short-period waves inside the bay (Pantoloan, Fig. 1b) and rapid arrival time onshore may indicate that tsunami might have been generated inside the bay by a relatively small source. This wave severely impacted on infrastructures and buildings at the front row near by the coastline (Fig. 2), however, the waves were not strong enough to penetrate further inland as the inundation line was found not more than 400 m landward. A main challenge for the source studies of the Sulawesi tsunami is the lack of enough sea level observations. Only two tide gauge records are available; i.e., Pantoloan and Mamuju. For inside the Palu Bay, where the extreme tsunami waves of 6-11 m were observed, only one sea level record is available (i.e., Pantoloan). The water waves from local phenomena such as submarine landslides usually travel short distances and decay rapidly (Fritz et al. 2004;McFall and Fritz, 2016;Heidarzadeh and Satake, 2017a,b); hence, even if we assume that a submarine landslide was involved, it is uncertain whether the single sea level record has captured such short-period and short-lasting waves or not.
A combination of spectral analysis, numerical simulations and field data may lead us to two different sources/mechanisms for the tsunamis recorded in Pantoloan and Mamuju. For the Pantoloan record, the inferred source dimension of 3-4 km ( Table 2) could be from a small submarine landslide; approximately with the same size as that of the 1998 PNG landslide tsunami. The larger source dimension of 32.5 km (Table 2) inferred from the Mamuju record is in agreement with the size of the initial co-seismic crustal displacement shown in Fig. 6c, d.

Conclusions
The 28 September 2018 Sulawesi tsunami was analyzed using sea level data and numerical modeling of tsunami. Main findings are: 1. The maximum tsunami wave height at two tide gauge stations of Pantoloan (within the Palu Bay) and Mamuju (outside the Bay) were 380 and 24 cm, respectively. The respective dominant wave periods were 3.6-4.4 and 10.0 min. The highenergy waves at Pantoloan lasted for 5.5 h while they lasted for at least 14 h in Mamuju. The two observed waveforms were significantly different with wave amplitude and period ratios of *16 and *3, respectively. 2. Based on the results of spectral analyses of the tsunamis recorded inside (i.e., Pantoloan) and outside (i.e., Mamuju) of the Palu Bay, the dimensions of the tsunami sources are estimated at 3.4-4.1 km and 32.5 km, for inside and outside of the Palu Bay, respectively. 3. Numerical simulations of the tsunami using a purely-tectonic source model based on USGS were capable of fairly well reproducing both observation tide gauge records in Pantoloan and Mamuju; except for arrival time in Mamuju. Our simulations (excluding runup calculations) resulted in maximum coastal amplitudes of *1.5 m which is far less than the observed runup of 6-11 m. After acquiring highresolution bathymetry/topography data and having a full understanding of runup field data, the model needs to be re-run. 4. Evidence may indicate the occurrences of submarine landslides following the Sulawesi event which potentially contributed to and intensified the main tectonic tsunami; such as the actual sea surface disturbances captured by video recordings from onboard an airplane and a boat within the Palu Bay. Submarine landslides could be responsible for extreme local runup heights; however, marine geological data need to be acquired from the Palu Bay to confirm such hypotheses. We identify the southern part of the Palu Bay, around the latitude of -0.82 o S, as the most likely location of a potential landslide based our backward tsunami ray tracing analysis.

Acknowledgements
Two tide gauge records in Pantoloan and Mamuju were provided by the Agency for Geo-spatial Information, Indonesia (BIG) (http://tides.big.go.id). We sincerely thank scientists in BIG including Andi Eka Sakya, Hasanuddin Z Abidin and Mohamaad Arief Syafi'i for their efforts to provide the data in a timely manner. The Latu Datu tide gauge data were downloaded from the Intergovernmental Oceanographic Commission website (http://www. iocsealevelmonitoring.org/). The earthquake source model used in this study is based on that of USGS (https://earthquake.usgs.gov/earthquakes/eventpage/ us1000h3p4/finite-fault). We are grateful to scientists in the USGS Earthquake Hazards Programs for making the source model publically available. We used the GMT software by Wessel and Smith (1998) in this study. We sincerely thank Alexander Rabinovich (Editor-in-Chief) for reading and commenting on the early version of this article and for producing Table 1. Results of spectral analyses were discussed with Alexander Rabinovich. We received help from Satoko Murotani (National Museum of Nature and Science, Tsukuba, Japan) for the analysis of vertical water displacement due to steep slopes. We sincerely thank two anonymous reviewers and Alexander Rabinovich (Editor-in-Chief) for their constructive review comments. MH is grateful to the Brunel University London for the funding provided through the Brunel Research Initiative and Enterprise Fund 2017/18 (BUL BRIEF).
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http:// creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.