Modelling of potential tsunami in Makassar strait, Indonesia: Submarine landslide-induced tsunami simulations

Tsunami modelling of potential landslide-induced tsunami in Makassar Strait is carried out to quantify possible damage to the nearby cities. Two numerical models are used to represent the wave generation and propagation by using NHWAVE and FUNWAVE models, respectively. The simulations consist of a series of scenarios based on distinct size of the landslide volume. Four landslides with volume 5, 8, 70, and 200 km 3 are used as tsunami sources in the initiation stage. The sources are evenly distributed in the Strait addressing different landslide location. Maximum wave heights of 1.5 m are found in the area between Palu and Bangkir from case 1 and around Talok from case 2 simulations. The empirical run-up calculation of 7.5 m is estimated at the land for the presented wave height. The value signiﬁcantly elevates the case 3 and 4 proportional to the volume values. The waves impact more than half of coastline with maximum value found in the Sulawesi side. Interestingly, wide and narrow shelf next to Kalimantan Island plays an important role in reducing the tsunami hazard level.


Introduction
A tsunami is a major threat to Indonesia archipelago due to its location surrounded by the active subduction zones. The tectonic stresses that are released from the subduction of Indo-Australia plate beneath the Sunda shelf are potentially generated tsunami 1 . For instance, one of the most devastating catastrophic occurred in December 2004, which triggered by M w = 9.3 Sumatra-Andaman earthquake. This tsunami generates 30 m run-up height in Aceh, Indonesia 2 and approximately 20 m in Thailand 3 .
Despite the majority of large tsunamis are produced by underwater earthquakes 4 , the latest tsunami event in Indonesia revealed the importance of studying other tsunami sources (such as subaerial, submarine landslide and earthquake-induced landslide). For instance, the flank collapse of Anak Krakatoa on 22 December 2018 generated tsunami and affected certain areas around Java and southern Sumatra, killed 437 people 5 . Field survey reported that this tsunami produced wave run-up more than 85 m at the northern coast of Rakata and 83 m at the southern coast of Sertung 6 . Rakata and Sertung islands located less than 5 km from Anak Krakatoa and facing directly its Southwest flank. A numerical study to investigate the effect of landslide volume in Anak Krakatoa tsunami has been done by 7 . By using an assumption that a rigid body with a volume of 0.15 km 3 sliding into water as a granular flow under gravity forces, destructive wave of 1.7 m height is observed in the Tanjung Lesung (approximately 50 km south east of Anak Krakatoa) after 28 minutes of propagation. This phenomena is described as the most damaging volcanically-resulted tsunami since the eruption of Krakatoa in 1883 8 . It's worth to note that the Anak Krakatoa tsunami shows that the tsunami induced by a subaerial or a submarine landslide can generate destructive tsunami waves 9,10 .
Generation mechanism of tsunami induced by a landslide is more diverse compared with the earthquake-tsunami, in which composed of the impulsive waves due to subaerial landslides falling into the water with high velocities [11][12][13] to huge submerged landslide 14 . Furthermore, a combination of earthquake triggers subaerial or submarine landslide-induced tsunami has elevated its threat to the coastal area. Particularly, submarine landslides can be triggered by a moderate magnitude of earthquake which take place in the continental slope 15 . The example of such mechanism occurs on 28 September 2018 tsunami in Palu, Sulawesi. An earthquake of 7.5 M w struck at the land of Sulawesi followed by a destructive and deadly tsunami producing measured run-up of 9.1 m at Benteng village. Ten large coastal sectors have reportedly collapsed into the sea after the earthquake.
The measured tsunami data collected around Palu Bay clearly show the contribution of secondary non-seismic local sources generating tsunami 16 . Similar finding has been reported in another research 17 . They concluded that among the arrived tsunami waves there are two initial waveforms most likely generated by a submarine landslide with approximately volume of 0.02-0.07 km 3 at the Southwest part of Palu 18 . As a result of 3.2 million m 3 of materials lost from the seabed causing a maximum decrease in the seabed elevation of 40 m. These results were supported by 19 which highlighted that the generation of Palu tsunami was a combination of sea floor changes as a result of earthquake followed by landslides. More importantly, tsunami heights from a combination of small to large submarine landslides along the Palu shores could reach up to 7.0 m 20 .
Due to its strongly varying initiation of motion and complex interplay between different triggers mechanism, the quantification of the tsunamigenesis from the submarine landslides are highly important. Using bathymetric data collected from MERAMEX (2004) and SINDBAD (2006) surveys, six submarine landslides in the eastern part of Sunda margin which span from central Java to Sumba Island have been identified 21 . The volume ranges from 1 km 3 in the Java basin to 20 km 3 in the Sumba trench. A series of simulations by using the aforementioned landslides unveiled run-up height of 7 m, 6 m, 3 m, and less than 2 m at Sumba, Lombok, Bali, and Java, respectively.
While major effort has been invested to study 22 and developing a tsunami early warning system, mitigation plans, and significant amount of works in tsunami modelling and risk assessments in the western part of Indonesia, less attention has been given to the eastern part of Indonesia. Numerical simulation of potential tsunami shows that the maximum tsunami height is more than 10 m locally and exceed 20 m next to the source 23 . This condition strongly indicates that the eastern part has similar order of tsunami threat compared with the western side of Indonesia. It should be noted that, more active regions strongly describe that the tsunami hazard and risk in the eastern side are alarmingly high.
Among the eastern area, the highest frequency of tsunami events for Indonesian archipelago takes place in Makassar Strait. A collaborative research has evaluated slope stability in the Makassar Strait as a result of a huge amount of sediment being transported from the Mahakam Delta 24 . The strait, in which is being known as the pathway for transferring water mass from Pacific to Indian Ocean can trigger slope failure due to massive sedimentation or erosion. Approximately 100-600 km 3 of alluvial sediment were observed on the west side of the Makassar Strait. This is a tectonically active region and a place where the sediment deposited. Previous tsunami study has revealed the necessity to account the secondary generation mechanism, submarine landslide 25 . In his research, submarine landslides are believed has strong contribution to add more energy of tsunami. Another work shows similar agreement indicated the importance role of landslide-induced tsunami in the Makassar Strait 26,27 .
The Makassar Strait separates western and eastern Indonesian archipelago. This area is classified as a marginal basin occupying the continental slope and rise regions within Sulawesi and Kalimantan Islands 25 . It's located at the intersection of four major plates, such as the Indo-Australian Plate in the southern part, Eurasian Plate lied in the west, Pacific Plate in the east, and Philippine Sea Plate in the north-west area. Developing a complex system of subduction, back-arc-thrusting, extension and major transform zones 24 . The seabed features of Makassar Strait distinctly describes two steep gradient zones as the southern and northern boundaries representing the Pastenoster and Palu-Koro transform fault. Historical data describe that seismic activities are mostly dominated by shallow depth earthquakes as a result of back-arc-opening zone in the strait and north-south Palu-Koro fault movements 1 . As expected, a previous study unveils lower magnitudes (5.5-6.6) of tsunamigenic earthquakes at the Pasternoster than at the another faults (6.3-7.7) 25 . He concluded that both empirical and numerical calculations of initial tsunami wave coming from Pasternoster sources are to small for destructing the coastline. Hence, a secondary mechanism such as submarine landslide is important to consider. In fact, Due to its narrow width, it is possible that far-field earthquake is the triggering mechanism for the landslide inside the strait.
A high rate of sediment transport results in under-consolidated accumulation of the Mahakam Delta is considered as the prominent contributor to the slope instability. It is fine grained sediment composed of high organic contents and gas charge, which significantly recede the sediment shear strength produced instability slope condition. Furthermore, interpretation of seismic, bathymetry and gravity core offshore the Mahakam Delta indicates high susceptible landslide. Following the aforementioned findings, the primary focus on this research would be on tsunami generation and propagation which account for submarine landslide. To the best authors knowledge, this study is the first landslide-induced tsunami simulation in Makassar strait.

Tsunami Simulation
The characteristics of tsunami generated by submarine landslide are mainly determined by the volume, initial acceleration, maximum velocity, and possible retrogressive behaviour of the landslide 28 . Tsunami generation in this research is simulated using a three-dimensional shock capturing Non-Hydrostatic Wave (NHWAVE) 29 . The model domain covers most part of the Makassar Strait (see Figure. 1) composed of uniform rectangular grids. Eights layers are imposed in the vertical direction. The horizontal grid resolution is 100 m yielding in total 277 grids in x-direction and 295 grids in y-direction. Four tsunami sources were distributed in the Strait with the total volume 4, 8, 70, and 200 km 3 following geographic, oceanographic, and seismic study 24 . These sources are modelled as a rigid body falls into the water and it excludes the deformation effects. The simulation time for each source is 30 minutes.
After tsunami characteristic time (t c ) is achieved in NHWAVE, simulation is continued by using less computationally demanding model, the Total Variation Diminishing (TVD) of the fully non-linear Boussinesq wave model (FUNWAVE) 31 . The FUNWAVE model uses larger domain than the NHWAVE domain, which fully presented sea region in the Strait. Uniform grid with 500 m resolution is used. The driving forces such as surface elevation, velocity component in x-direction, and velocity component in y-direction in the FUNWAVE model are obtained from NHWAVE solutions after t c is fulfilled. This is a common procedure in modelling landslide-induced tsunami due to the multiple spatial and temporal scales being involved 32 . Bathymetry data with 6" resolution used inside the both domains were obtained from National Hydrographic Chart (BATNAS). In addition, detailed information about the landslides can be seen at Table 1. Analysing of landslide-induced tsunami simulation outputs presents deep insight into tsunami characteristic at the Palu-Koro fault. During the initiation stage, a dipole wave with circular shape formed over the failure area ( Figure 2a). It spreads radially and progressively deforms more circular wave. Positive waves propagate to the South while negative waves move to the north. Wave heights were in the range -1 to 1 m outside the centre mass and peak values reached to about 5 m. The tsunami waves maintain its radial shape during the penetration to and out of the Strait as shown in Figure 2b.

Results
The wave trains propagation to the south occurs at lower speed, greater amplitude and shorter wavelength (Figure 2c). These waves impacted the Sulawesi Island 20 minutes after landslide initiation. The first reporting positive wave was at Bangkir synthetic wave gauge with amplitude of 0.7 m, in which followed by -0.8 m wave trough. The Figure 3b1 clearly describes that subsequent waves arrived periodically at 10-12 minutes interval with significantly lower amplitudes. After approximately 120 minutes since it's generation, the tsunami waves were observed at Palu synthetic gauge. Wave amplitude of 0.15 m and 10-12 minutes wave period were recorded in this site (Figure 3b2). The tsunami waves required more time to reach the Kalimantan Island compared to the Sulawesi Island. The first-arrived tsunami wave at the Talok synthetic wave gauge was the minimum surface elevation of -0.2 m followed by the maximum free surface elevation of 0.15 m as shown in Figure 3b4.
The maximum surface elevation plot (see Figure 3a) shows the maximum value of surface elevation in each grid during whole simulation time. For the sake of clarity, the maximum surface elevations less than 0.3 m were removed from the plot. The values were in the range of 0.8 to about more than 5 m in the simulation domain. The maximum value of 1.8 m was observed in the area between Bangkir and Palu with the tsunami time arrival less than 30 minutes. In the southern part of Palu, the maximum surface elevations of 0.3 -0.8 m were recorded with the arrival time approximately 50 minutes. Analysis from the maximum surface elevation plot explains that there is no coastal region in the Kalimantan Island observed surface elevation more than 0.3 m (see Figure 3b4). Its worth to note that, the plot not only unveils the wave coming from direct tsunami but also the tsunami reflections, diffractions and refractions. The initiation stage of tsunami wave generation is shown in the Figure 4a. Similar with the case 1, a dipoles observed at the initiation stage. Less speed, greater amplitude, and shorter wavelength tsunami waves penetrate to the south where deposited sediment occurred. Contrastingly, a tsunami trough wave moved to the north with greater speed before being reflected back to the Strait. The wave crest and trough initial value after the landslide occurred are 3 and -6 m at the centre mass. After 25 minutes from tsunami generation, the wave trains started to reach shallower depths while still able to maintain its circular shape (Figure 4b). The first-recorded tsunami wave at Bontang station was -0.14 m of trough wave about 50 minutes after initiation (see Figure 5b3). These waves was then followed by 0.13 m high wave 10 minutes later. The Palu station observed a 0.4 m high wave 90 minutes after the landslide as shown in Figure 5b2. During this time, the tsunami wave has propagated far away to the south of the Strait (Figure 4c).
The maximum elevations plot distinctively shows higher value at the southern part of Palu coastline (Figure 5a). Interestingly, the maximum values are in the range of 0.6-0.8 m, in which two times higher than other coastline areas (less than 0.3 m). We believe that bathymetry condition around Palu Bay was the key point in this condition. The average value of the tsunami time arrival in this region is 90 minutes.

Case 3: Submarine Landslide next to Mahakam Delta
The scenario is a failure triggered by massive sedimentation in Mahakam Delta. The total volume of the landslide is 70 km 3 following previous study of 24 . Similar landslide parameters with the previous cases are applied.
The landslide initially produced 25.2 m of high and -32.0 m trough waves that propagate perpendicular to the shore as highlighted in Figure 6a. After 25 minutes simulation, the wave train started to reach shallower depths. Maximum wave with 4.89 m high propagated to the east while negative wave of 3.2 m moved to the west with period of 10 minutes (see Figure 6b). The first-observed tsunami wave was at Bontang station with -2.0 m wave trough 50 minutes after the landslide, followed by 3.0 m high wave approximately 10 minutes later as shown in Figure 7b3. The amplitude of this wave is much higher than the previous two cases. Then, a 1.5 m positive wave was recorded at Palu station 70 minutes after generation time (Figure 7b2), while at Talok station 0.2 m high wave was observed approximately 50 minutes later (Figure 7b4). Propagations to the north still continue even after 100 minutes from the initiation time.
To address the maximum wave elevation close to the coastline area, the maximum value of surface elevation is analysed. The Figure 7a

Case 4: Worst case scenario of submarine landslide at Mahakam Strait
In this scenario, a total of 200 km 3 sediments is used as the driving force for tsunami generation. The number of sediment volume is derived from comprehensive study by considering sediment deposit coming from Mahakam Delta 24 . As a result, the landslide is categorised as the failure due to the massive sedimentation and moved as a translational slide. This slump occurred at 1500 m water depths, with 0.25 m/s 2 initial acceleration, and 67.5 m/s terminal velocity. The value of terminal velocity is equal with the Storega landslide simulation and it's assumed as the possible value of maximum velocity. The time characteristic produced by this scenario is 270.0 m/s. Figure 8a illustrates the initial tsunami wave generation mechanism adjacent to the Mahakam Delta producing approximately 40 m surface elevation at the centre mass and 0.5 elevation at the outer side. The wave trains moved in the west-east direction and took form of circular waves (Figure 8b). While propagate to the West, the waves reduced its speed and wavelength, but increased its magnitude as the bathymetry changed from deep to shallow depth. Contrastingly, as the wave move to the East, the speed and wavelength elevated but decreased its magnitude. Amplified surface elevation is observed after more than 1 hour simulation as the wave commenced to reach the Shelf of Samarinda and Southern part of Palu coastline (see Figure 8d and Figure 9b4). The wave trough with magnitude less than 0.05 m was observed at Bontang station after 110 minutes which followed by 0.18 m high wave 10 minutes later as shown in Figure 8b. Higher wave amplitudes were recorded at Samarinda and Palu station after 120 minutes. The first arrived waves at Samarinda station was -0.5 m trough wave followed by 1.5 m high wave 10 minutes later. These waves clearly showed typical of tsunami waveform. Figure 9b2 describes that high wave with 2.0 m amplitude is found at Palu station. The second observed wave at this station is -3.8 m wave about 10 minutes after the first wave. The following wave trains at Palu station are still high, even after the second wave with 1.0 m wave magnitude as indication of amplification process due to wave reflection, undulation, and refraction.
The maximum surface elevation showed disparity values in the range of 3.0-9.3 m, particularly in the side of Sulawesi Island (southern part of Sulawesi) as clearly depicted in Figure 9a. These sites are directly hit by the waves during simulation. Interestingly, there are no areas with surface elevation maximums exceed 0.3 m in the side of Kalimantan Island despite some waves propagated directly to the shore. Overall, the patterns are similar with the case 3 where the amplification at Sulawesi Island exist.

Discussion
An initial standpoint of submarine landslide as secondary tsunamigenic mechanism in Makassar Strait has been previously reported 25 . In their report, seismic profile data of the Southern Makassar Strait described the condition of an older trench filled with slumping materials. Unfortunately, there is no further study to unravel the tsunamigenic sources in the Strait until recently geological study of Makassar strait is published 24 . By using comprehensive data comprised of bathymetric, oceanographic, and seismic, their research concluded that significant amount of sediments has been deposited along the shelf next to the Mahakam Deta.
The analysis of tsunami hazard due to the landslide failure in this study is a scenario-based tsunami hazard assessment (STBHA) since there is no primary data of landslide in this region. Such procedure has shown satisfying result in some instances, if the imperfection and uncertainty in the assumed scenario is properly considered. The assumptions in this research are there is only single landslide for each scenario (i.e. there is no subsequent slump during simulation) and only accounts for four landslides. Another critical point is that no direct run-up calculation originating from the model because of poor resolution of the grid. Horizontal grid size at least 50 m is required to accurately resolved computation of tsunami inundation and run-up 34 . Hence, the run-up calculation in this study used empirical formula 35 , in which significantly influenced by slope condition at nearshore and wave height at continental basin.
In order to understand the threat from tsunami to the surrounding coastline, it is important to provide a tsunami wave height threshold. Experienced from previous tsunamis 2018 in Palu and Krakatoa, extreme damage to houses and vehicles may occur due to strong water mass flux and currents in the nearshore. The minimum tsunami wave amplitude or run-up at 1.5-2.0 m has ability to caused severe damage 36 . Another study chose tsunami run-up of 1.5 m for potential coastal damage and 3 m for having significant consequence to the coastal 37 .
The simulations carried out in this study indicate potential severe damage to the local coastal area in Makassar strait. Impacts are varied depend on the generating mechanism, landslide characteristics and regional bathymetry. Simulations with landslide volume less than 10 km 3 (case 1 and case 2) reveal few sites with wave more than 1.5 m high. These areas are the Southern side of Bangkir (case 1) and the bottom part of Kalimantan Peninsula (case 2). By assuming the slope in that region is 10 o and 5 m of water depth (d), the calculated run-up from 35 's formula is 7.5 m which is believed has enough energy to produce severe damage at the land. Notice that, these are resulted from mid-size landslide and even less than 1/10 volume of the worst scenario. The third and worst case simulations will significantly elevate the hazard of tsunami to the both islands.
More than half of coastline in the North Makassar Strait affected by the waves according to the simulation results. The maximum tsunami waves are in the range of 1.5-3.0 m in the Kalimantan side and higher values are recorded to about 3.0-7.8 m in the Sulawesi Island. From the simulations, an interesting fact is observed where the tsunami waves propagated to the East (toward Sulawesi) gained more energy than the waves moved to the West. As a result, the simulated wave heights in the Sulawesi side are frequently larger than at Kalimantan side. This result is consistent with the simulated wave height from 25 which caused by bore formation of the leading wave on the Kalimantan coast due to the large, relatively shallow shelf area.
The landslide tsunami presents more challenges compare with the tsunamis generated by earthquake. For instance, the subduction faults are generally well known and has been studied by broad communities. By comparison, landslides generated tsunami are nature local features and very difficult to be uncovered except with collective research. Since Palu tsunami 2018, the government has deployed many sensors around the archipelago, but still there is no particular technologies to deal with the landslide-induced tsunami. As a result, a number landslide-tsunami simulations is required to provide more understanding and alert to the society.

Conclusions
Tsunami generation and propagation simulations are conducted in order to analyse the potential tsunami hazard in Makassar Strait caused by a submarine landslide. The Strait is defined as a highly tectonic active region and its location where a significant amount of sediment being deposited. The sea bed features of this area clearly show two alignments of the continental shelf in the side of Kalimantan and Sulawesi Islands. Based on the aforementioned condition, both deposition or erosion has possibilities to trigger slope failure and due to the narrow width, far-field earthquake can also destabilise the slope.
The simulation results clearly show typical waveforms of landslide-tsunami identified by a high wave propagates toward the slump direction. Different wave heights are produced from the simulation, which depend on the volume size demonstrating various damage levels to the shore. A landslide with volume less than 10 km 3 is calculated and produced 0.8 to 1.5 m wave height at 5 m depths. Estimated run-up of 7.5 m is obtained for 1.5 m wave heights. Case 3 and case 4 produced maximum tsunami wave in the range of 1.5-3.0 m in the Kalimantan side and 3.0-7.8 m in the Sulawesi side. As the result, higher run-up and stronger tsunami waves are predicted from the two worst scenarios.
Due to the lack of landslide parameters in Makassar Strait, a scenario-based study is used which underlie on the assumption that only single source of the landslide exists for each scenario. There is no subsequent slope failure in the simulation. The assumption might lead to under-predicted results. Therefore, it is important to conduct a comprehensive survey as a way to provide primary data. Case Figure 1. Study location comprises of NHWAVE and FUNWAVE model domains. Bathymetry and topography data were obtained from The Geospatial Information Agency of Indonesia (BIG). Landslide locations were defined based on geographic, oceanographic, and seismic study 24 (grey diamonds). Earthquake data were retrieved from Global CMT https://www.globalcmt.org/(red circles).