Solving the puzzle of the 1996 Biak, Indonesia tsunami

On February 17, 1996, an earthquake ( M w = 8.2) occurred northeast of Biak Island, Indonesia, and generated a tsunami. Interestingly, the tsunami runup on the southwest side of Biak Island, which did not face the epicenter, was higher than that on the side of the island facing the epicenter. It was earlier suggested that the earthquake triggered submarine landslides. However, this hypothesis remains unexplored and unconfirmed. In this study, tsunami arrival times obtained from eyewitness accounts were used to perform backward tsunami travel time modeling. Based on the backward tsunami travel time results and multibeam bathymetry, five submarine landslide candidates were identified: three located to the southwest of Biak Island and two to the south of the island. Our analysis indicates that tsunamis from a large submarine landslide off the southwest coast of Biak Island (SL 2) and a small sub-marine landslide off the south side of Biak Island (SL 4) contributed to the 1996 Biak tsunami event and came earlier than the main tsunami from the earthquake fault. Tsunami sources were earlier modeled only by fault parameters without considering submarine landslide sources. As a result, the models could not explain the observed runup and eyewitness tsunami arrival times for the southwest and south of Biak Island. To address this problem, we propose an approach that combines a submarine landslide model with a modified version of a previously proposed fault model. Our combined model explains the observed runup and eyewitness tsunami arrival times well (Aida’s index values of the model are K = 1.00 and κ = 1.44).


Tectonic setting surrounding Biak Island
Biak Island is located in eastern Indonesia to the north of Yapen Island (Fig. 1).This region encompasses several faults that generate seismic activity.The Sorendidori fault, which runs directly beneath the island, divides Biak Island into two parts: the Supiori district in the north and the Biak Numfor district in the south.It is a strike-slip  (Ekström et al. 2012) fault with a movement of 0.5 mm per year (Pusat Studi Gempa Nasional 2017).The Biak Fault Zone (Fig. 1) is present in the southwestern part of Biak Island.This zone comprises a series of parallel strike-slip faults which run along the southwest coast of the island (Gold 2018).The New Guinea Trench (NGT) is located to the north of Biak and Papua Islands.This trench was formed by the collision of the Australian and Pacific Plates.In the west and south of Biak Island, there are the Manokwari Trench and the Yapen Fault.According to earthquake data from the International Seismological Centre (2023) and the Indonesian Agency for Meteorology, Climatology, and Geophysics (2023), the seismic activity around Biak Island was not very high in the period 1975-2022.Generally, large earthquakes M > 5 occurred to the east of Biak Island and were associated with the NGT (Fig. 1).

Graphical abstract
Based on historical data from the National Geophysical Data Center/World Data Service (2023), four tsunamis occurred near Biak Island in 1914Island in , 1979Island in , 1996Island in , and 2009 (Fig. 1).The 1914 and 1979 tsunamis were associated with earthquakes of magnitudes M w = 8.1 and M w = 7.9, respectively.The 1979 event was attributed to the Yapen Fault.However, it is unclear whether the Yapen Fault or the NGT caused the 1914 event (Okal 1999).The 1996 tsunami was caused by the NGT (Okal 1999).The most recent tsunami occurred in 2009 and was generated by an earthquake of magnitude M w = 7.7 centered on the Manokwari Trench (Fujii et al. 2011).

The 1996 Biak tsunami
On February 17, 1996, at 05:59 UTC, an earthquake of magnitude M w = 8.2 occurred northeast of Biak Island and generated a tsunami (Fig. 1).Interestingly, this was a thrust-fault earthquake which occurred in the NGT segment, which has no previous historical record of large earthquakes (Okal 1999).The tsunami did not only hit Biak Island but also reached Taiwan, Japan, Canada, and the USA with maximum heights of 0.28 m, 1.04 m, 0.18 m, and 0.35 m, respectively (Imamura et al. 1996;Li et al. 2006;Stephenson et al. 2007;Lau et al. 2010).One month after the event, the International Tsunami Survey Team, comprising tsunami researchers from various countries such as Japan, Indonesia, and the USA, conducted a tsunami survey on Biak Island (Imamura et al. 1996(Imamura et al. , 1997)).They found that the tsunami hit the entire coast of Biak Island and had an average height of 2 m.Surprisingly, a tsunami with a height of 7.7 m hit Mardori in the southwestern part of Biak Island.This location did not face the epicenter of the earthquake, but the tsunami runup here was higher than that on the northeastern coast, where the maximum runup was 5.4 m at Korem (Fig. 2).This event left 110 people dead, 51 missing, and 100 injured; about 10,000 people lost their homes.The damage was the most severe in a small village in a narrow bay in Korem, where the tsunami wave struck and swept away all houses at a distance of up to 800 m from the shore (Imamura et al. 1996(Imamura et al. , 1997)).
This tsunami was simulated by Imamura et al. (1997) using a fault-only model.They used the moment tensor solution from the Harvard University Centroid Moment Tensor to set the fault parameters.The results of this simulation were mostly in agreement with the observation data, except for Mardori.For Mardori, the calculated runup was only about 1 m, but the observation data indicated a height of 7.7 m.Therefore, they suspected that the M w = 8.2 earthquake northeast of Biak Island trig- gered a local submarine landslide (SL) on the southwest side of the island.Matsutomi et al. (2001) also suggested the possibility of an SL.They interviewed the local populace about the arrival time of the tsunami and concluded that the tsunami wave hit the south and southwest coastal areas earlier than the northeastern coast possibly because of the occurrence of many SLs near the south and southwest coasts.Another case of superposition of waves coming from different sides at the back of an island was reported by Heidarzadeh and Satake (2013) with reference to the Mediterranean tsunami that occurred in May 2003.
Several tsunami surveys, damage assessments, and runup simulations were conducted to study the Biak Island tsunami (Imamura et al. 1996(Imamura et al. , 1997;;Matsutomi et al. 2001).However, a satisfactory explanation for the unexpectedly high runup in the southwest coastal area is still lacking.Until now, the SL which contributed to the 1996 Biak tsunami remains only an assumption.Therefore, in the present study, we attempted to prove this hypothesis.

Bathymetry and topography
Three nested grid layers with resolutions of 300 m (Layer 1), 100 m (Layer 2), and 50 m (Layers 3a and 3b) were used in this study (Fig. 2).A combination of bathymetry and topography data (BATNAS; Batimetri Nasional) was used in Layer 1.These data can be downloaded from the website of the Indonesian Agency for Geospatial (2024).High-resolution multibeam data were procured from TGS (UK); these data cover the southwest sea area but not the shallow coastal area.To obtain bathymetry information not provided by the multibeam or terrestrial data, we used nautical chart numbers 224 and 225 procured from the Indonesian Navy (2023).The Digital Elevation Model Nasional (DEMNAS) provided by BIG is the primary source of topography data used in this study; these data can be freely downloaded.The nautical chart, multibeam, and DEMNAS data were combined for Layer 2 and Layers 3a and 3b (Fig. 2).

Tsunami runup
The runup height (Imamura et al. 1996) was determined based on a field survey (Fig. 2).The runup data covered the entire area of Biak Island, except the northern part.In the northeastern coastal area, facing the epicenter, the highest runup height of 5.4 m was recorded at Korem.The tsunami runup in the northeastern coastal area had almost the same height and then it gradually decreased in the northern part of the island.Meanwhile, the tsunami runup increased significantly in the southwest coastal area, which did not face the earthquake epicenter and had an inverted V-shaped coastline.The highest tsunami runup recorded in this area was 7.7 m at Mardori.
To investigate the possibility of an SL, we applied the same technique as Okal et al. (2002).We performed fitting using a weighted least square with order 2 for the northeast, southwest, and southern coastal areas using the "trend1d" command in Generic Mapping Tools 6.3 (GMT; Wessel et al. 2019), together with manual smoothing (Fig. 2).The results are compared with those of Okal et al. (2002) in Fig. 2. The tsunami runup distribution for the northeast and southern coastal areas has the same pattern as that of the 1992 Nicaragua tsunami, while the southwest coastal area has a profile similar to that of the 1998 Papua New Guinea tsunami.This result suggests that the northeast and southern coasts were affected by a tsunami generated by the earthquake itself, while the tsunami caused by an SL hit the southwest coast.

Eyewitness tsunami arrival time
The tsunami arrival time is defined as the initial instance when the water level is changed.For eyewitnesses, the arrival time is defined as the moment when they first observed a noticeable change in sea level after the primary earthquake shock.Despite potential variations in these arrival times, we consider them under the same definition.This comparison is solely used to confirm the approximate location of the tsunami source area.Matsutomi et al. (2001) estimated the tsunami arrival time for the Biak Island by interviewing the local inhabitants (Fig. 2).Their questionnaire included questions on characteristics of the incoming waves, such as the period, number of waves, and highest wave.For the area facing the epicenter, the eyewitness tsunami arrival time was 10-15 min after the mainshock.For the southwest and southern sides of the island, the eyewitness arrival time was ~ 5 min after the mainshock, and the shortest eyewitness arrival time was 1 min after the mainshock at Mardori.The differences in these arrival times seem to indicate another source of the tsunami.

Backward tsunami travel time
Many researchers have used the backward tsunami travel times (TTT) method in identifying the source of tsunamis.Heidarzadeh and Satake (2014) and Heidarzadeh et al. (2023) applied this method to identify the sources of tsunamis resulting from inland earthquakes in Pakistan and Turkey, respectively.In this study, we used this method to locate any possible SLs in the southwest and southern coastal areas of Biak Island.The eyewitness arrival time of the tsunami, as reported by the local people (Fig. 2), was used to plot the backward TTT lines.Since accurate location information could not be obtained from the eyewitnesses, we assumed that the observation points were near the coordinates of the tsunami runup survey (Fig. 2).For travel time estimation, we used the TTT software package (TTT SDK v 4.0.1)provided by the International Tsunami Information Center, assuming long-wave approximation based on Huygens' principle (Wessel 2009).

Tectonic source
Since an earthquake was the main cause of this event, some researchers tried estimating the fault parameters.Imamura et al. (1997) and Tanioka and Okada (1997) calculated the fault parameters using the scaling law.Henry and Das (2002) obtained the slip distribution by inverting both P and S waves from seismographs, while Arimuko and Sianipar (2021) only used P or S waves.Unlike other fault models, the model proposed by Henry and Das (2002) was not completely rectangle (Fig. 3a).
In this study, the abovementioned fault models were used in tsunami simulations to determine whether they can explain the tsunami runup distribution on all sides of the island.Since Arimuko and Sianipar (2021) did not provide detailed subfault parameters, for simplicity, we assumed that their model represents a single fault.We also set a depth of 10 km for the fault reported by Tanioka and Okada (1997), because they did not specify the depth value.All fault models were simulated instantaneously, except for the Henry and Das (2002) model.The Henry and Das (2002) model ruptured from S2-S4-S3-S1.The details of all fault parameters used in this study are given in Fig. 3a and Table 1.

Tsunami simulation
The tsunami that reached the southwest and southern coastal area is assumed to have mainly resulted from an SL (Fig. 2).A sliding mass causes the initial sea surface wave to be approximated as a double three-dimensional Gaussian function (Watts et al. 2005).To estimate the wave height, we adopted a method based on a slumptype mechanism (Watts et al. 2005).This method was earlier used to explain the SL during the 1771 Ryukyu tsunami (Okamura et al. 2018), the 1992 Flores tsunami (Pranantyo and Cummins 2019), and some SL events in eastern Indonesia (Pranantyo et al. 2021).
To calculate the tsunami propagation from the source to the coast, the nonlinear long-wave equations in Cartesian coordinates were solved using the leapfrog finitedifference scheme (Goto and Ogawa 1997).We used the Boussinesq-type model with a dispersive effect for Layer 3, since an SL creates a tsunami that affects only a small area and generates shorter wavelength waves that are more dispersive (Løvholt et al. 2008;Fuhrman and Madsen 2009;Shi et al. 2012;Tappin et al. 2014).
A finer grid is not useful in the simulation, because nearly half of the observed runup coordinates were incorrect, because the data corresponded to the shallow sea area rather than the land.From 1990 to 2000, for military reasons, the US military intentionally introduced errors into the data for civilian Global Positioning System (GPS), a process known as Selective Availability.This practice reduced the accuracy to approximately 100 m (GPS 2000).To avoid this problem, it became necessary to simplify the comparison between the observed data and the numerical results.We selected a maximum resolution of 50-100 m for the nonlinear simulation (Layer 2 = 100 m and Layer 3 = 50 m in Fig. 2), as this resolution is sufficient to reproduce historical tsunami heights (Grilli et al. 2019); then, we selected the calculated runup at the nearest settlement or land.We allocated 1 h for tsunami propagation and a computational time step of 0.05 s to satisfy the stability condition.

Result validation
Aida's K and κ indices (1978) are widely used in validat- ing tsunami simulations (Koshimura et al. 2009;Nakamura 2009;MacInnes et al. 2013;Gusman et al. 2014).This technique is crucial for understanding the statistical significance of model accuracy and uncertainty.In this technique, the geometric mean, K, and the geometric standard deviation,κ , of the tsunami wave heights are used to evaluate the results.K is the mean value of K i (ratio of observation to simulation results determined at different points i), while κ is the mean standard devia- tion between K i and K.Because these indices are based on ratios (multipliers), it is necessary to use the geometric mean and geometric standard deviation instead of the arithmetic mean and standard deviation; however, the statistical understanding remains the same.By examining the geometric mean and geometric standard deviation of the ratio of the calculated and observed values, we can determine how many times, on average, that the model results differed from the observed values and the extent of the variation.Further, for the Aida (1978) score, the Japan Society of Civil Engineers (2002) set the threshold for a valid tsunami source at 0.95 < K < 1.05 and κ < 1.45.Thus, the parameters K and κ are defined as follows: where n indicates the number of data points, and x i and y i are the observation and simulation data at point i, respectively.

Geology of the Biak basin
The multibeam data revealed a basin between Biak Island and Yapen Island (Fig. 1 and Additional file 1).The predominant deposit of this basin is deep-sea clastic, and the basin has relatively low sediment contents.Carbonates form on the margins of this basin, and volcanic deposits are also seen locally and on neighboring islands (Bertoni and Álvarez 2012).A freshwater lens is thought to (1) , be the driving force behind the underwater spring-sapping process.Offshore Biak Basin, this lens is more than 55 km wide and about 250 m thick.Multibeam bathymetry revealed that spring-sapping by this freshwater lens caused several SLs in the Biak Basin (Gold 2018).In addition, the margin of the basin has a slope of around 5°-15°, which is enough to trigger an SL (Hance and Wright 2003;Bertoni and Álvarez 2012).
We analyzed and delineated the seafloor topography to identify the contour profile related to the SL, and we found many SL traces (Additional file 1).From the shallow to the deep seafloor, there is a canyon/gully where an erosional valley is formed by the small-scale collapse of the upper part of the slope.The Biak Fault Zone is also clearly visible (yellow dashed line in Additional file 1).Earthquakes with the strike-slip mechanism from the Biak Fault Zone were responsible for the occurrence of the SLs (Bertoni and Álvarez 2012;Gold 2018).

Tsunami simulation using a fault model
Two-layer nested grids (Layers 1 and 2) were used to simulate tsunami propagation with the fault-only model.This simulation was performed in determining whether all observed runups and eyewitness arrival time, especially the unexpected maximum runup of 7.7 m and 1 min eyewitness arrival time, can be explained by the fault models.
For all fault-only models (Table 1), the tsunami takes ~ 40 min and 25 min to reach the southwest and southern coasts of Biak Island (Additional file 2), respectively.This result is inconsistent with the eyewitness accounts according to which the tsunami took only about 5 min to reach the southwest and southern coastal areas (Fig. 2), indicating that this tsunami was not directly caused by the earthquake that occurred to the northeast of Biak Island.
All fault-only models can generally explain the distribution of the tsunami runup on the northeast and southern coasts (Fig. 3b, d).However, none of them can explain the 7.7 m tsunami runup on the southwest coast (Fig. 3c).Thus, these results strengthen the hypothesis of the existence of an SL (Imamura et al. 1997;Matsutomi et al. 2001).

Backward tsunami travel time results and submarine landslide locations
The runup distribution in Fig. 3c shows that the tsunami in the southwest area was caused by an SL triggered by the shaking during the mainshock or one of the aftershocks.The shortest eyewitness tsunami arrival time was only 1 min, and the ISC catalog showed that the closest aftershock with m b = 6.3 occurred about 1.5 min after the mainshock.We assume that the SL was triggered by the seismic wave of the mainshock and occurred immediately after that.
The SL was probably located near Mardori (Fig. 2).Bathymetry data with 50 m resolution in Layer 2 were used for the backward TTT calculation (Fig. 4).We did not use the data for Rayori and Owi, because in our calculations, the 5 min tsunami wavefronts were trapped in the narrow strait offshore of these points possibly due to the shallow depth and inaccurate bathymetry.No eyewitness information on tsunami arrival time is available for Wardo.Hence, this location could not be used for the backward TTT estimation.However, the sea receded first at this location (Matsutomi et al. 2001), and this information was useful to validate the SL simulation later.
Because of the uncertainties in the eyewitness arrival times and the locations of the eyewitnesses of the tsunami, the location of the SL could not be located accurately.Since the data were limited, we roughly estimated the possible SL locations using the backward TTT, and we add + 1 min buffer time to strengthen the consistency of the result for each eyewitness observation point (Fig. 4).Unlike the case of + 1 min buffer time, where the possible buffer area (red dashed circle) can be obtained from three backward TTT lines, we did not plot the − 1 min buffer time, because only two backward TTT lines intersected and the − 1 min buffer time lines also involve many uncertainties.
To determine how many SLs contributed to this tsunami event, we compared the distribution of tsunami runup in the southwestern coastal area (Fig. 2) with the data for the 2018 Palu, Indonesia, tsunami (Fig. 1 in Nakata et al. (2020)).In the case of the 2018 Palu tsunami, most SLs occurred near the east coast of Palu Bay and produced many peaks in the runup distribution.However, in the case of the 1996 Biak tsunami, in the runup distribution has only one peak near Mardori (Fig. 2).This suggests that only one SL contributed to this event.In addition, the shorter eyewitness arrival time in the southern part of Biak Island indicates another possible SL (Figs. 2 and 4).The backward TTT results indicate that this SL was located southeast of Samber.From the runup distribution pattern of the southern part of the island and the highest runup around Sorido, we predicted that the SL was small and faced around Sorido.
The detailed bathymetry cross sections are shown in Fig. 5. Inside the possible location of the SL and + 1 min buffer of backward TTT (solid and dashed green circle line in Fig. 5a), we identified three SL candidates on the southwestern coast (SL 1, SL 2, and SL 3).We did not choose the SL inside the canyon/gully, because it seemed small and could not trigger a high runup (Fig. 5a); we also excluded the SL to the right of the canyon/gully, because it was too far from Mardori.For the southern part, we identified two SL candidates (SL 4 and SL 5).Although the multibeam bathymetry did not cover the possible area (circle made of green solid line in Fig. 5b), we speculated that SL 4 was located at that point, because our assumption was also supported by the presence of a small topographic collapse moving southward to the south of the possible area.This analysis indicates that Sorido was located in the backside area of SL 4. Another possible SL source in the southern part was SL 5.This location was discovered from the + 1 min buffer time of the possible area (dashed green circle in Fig. 5b).The slight disturbance in the bathymetry contour in the purple rectangle and oval in Fig. 5b indicate the source and deposit of SL 5, respectively.

Submarine landslide parameter and simulation
The SL parameters of SL 1, SL 2, SL 3, and SL 5 were according to the geometry of the traces, because these SLs are clearly visible from the bathymetry data, but the parameters of SL 4 are hypothetical, because its location was not inferred from the multibeam bathymetry data (Fig. 5b and Table 2).For SL 4, if we extend the SL deposition trace to the northeast, the source geometry may not be larger than the possible area (green circle in Fig. 5b).In addition, since the SL 4 deposition site was not clearly visible, we assumed that it was transported to the SL 5 deposition site over a distance of 9500 m (Table 2).
Before attempting to identify which SL source in the southwestern area contributed to the event, we first compared Layers 2 and 3 to check whether these layers differ greatly in the calculated runup, since these layers have different resolutions.We use the fault-only models (Table 1) for the simulation, because these models are sufficient to generate a tsunami that hit the entire coast of Biak Island.The results showed that there is no large difference between Layers 3 and 2 in terms of the calculated  tsunami runups (Fig. 6).Then we calculated the root mean square error (RMSE) for all fault models for Layer 3a-Layer 2 is ± 0.11 and that for Layer 3b-Layer 2 is ± 0.16, with an overall average RMSE of ± 0.15.We simulated all southwest SL sources using the three layers (individually) to check the effect on the other side of the island before calculating the K and κ indices for the valid SL source.The initial sea surface elevations were determined to be about 0.8 m, 4.9 m, and 0.04 m for SL 1, SL 2, and SL 3, respectively.Our investigation showed that only small tsunami waves of SL 2 reached the southern coast and that no SLs hit the northeastern coast (Fig. 7a-c).This result is in agreement with the characteristics of SL/underwater slumps, which produce the highest runup and have little effect on the far field (Harbitz et al. 2006(Harbitz et al. , 2014;;Heidarzadeh et al. 2014).Then, we compared the calculated runup between SL 2 and all fault-only models on all sides of Biak Island.Most faultonly models show better calculated runup distributions than the SL 2 for both the northeastern and southern parts (Figs. 3b,d,7c), but for the southwestern part, SL 2 is better than all the fault-only models (Figs. 3c,7b).Hence, for calculating the K and κ indices for the south- west SL sources, the observed runups in the south and northeastern coasts can be ignored.The calculated runup from the three SLs at the southwest part are still underestimated (Fig. 7b) because of the limited and inaccurate observation data.However, the computational accuracy of our model results was significantly better than to the results obtained considering the fault-only models (Figs.3c and 7b).We compared the calculated tsunami arrival time with the tsunami characteristics reported by eyewitnesses, such as the eyewitness arrival time and first wave impulse (Table 3).We obtained good agreement for the first wave impulse for SL 2, but not for SL 1 and SL 3, which produced different results.Even for SL 3, the tsunami did not reach the eyewitness observation point in our simulation.Besides, with the limited data, it is difficult to reproduce the calculated tsunami arrival time according to historical records, especially the 1 min eyewitness arrival time.Moreover, we calculated the K and κ for the SL sources in Layer 3a to be 10.22 and 2.07, respectively, for SL 1 and 1.00 and 1.37, respectively, for SL 2. The K and κ results of SL 1 and SL 2 are contrasting, because the size of the SL 1 was really small compared to SL 2, and hence, SL1 could not trigger higher calculated runups (Fig. 7b).Thus, from this analysis, we believe that SL 2 contributed to the 1996 Biak tsunami event.
The calculations of K and κ indices for SL 4 and SL 5 were not essential, because these SLs produced small tsunami waves and did not significantly impact the land (Fig. 7a-c).Since these SL waves were very small on the other sides, we only focused on identifying the SL source which contributed to the tsunami which affected the southern coast.From Fig. 7c and Table 3, we concluded that SL 4 was the main source at the southern coast, because both the calculated runups and the calculated tsunami arrival times were better described by SL 4 than by SL 5.Although SL 4 contributed to the tsunami that affected the southern coast, the observed runups in this area were well described by the fault-only models (Fig. 3d) and showed a pattern similar to the northeast observed runup (Fig. 2).We assumed that the eyewitnesses might have seen the tsunami from SL 4 at that time, but the collected runups in the shoreline were due to the tsunami from the main earthquake.

Tsunami simulation using submarine landslides and fault model
To create a model suitable for the observed runups on all sides of Biak Island, we perform a tsunami simulation by combining the fault models (Table 1), SL 2, and SL 4 (Table 2) at the same time (Fig 7d-f ) and then estimated K and κ for each model combination.The best result was obtained by combining the fault model of Imamura et al. (1997) and SL sources, because it produced the smallest κ value and the smallest geometric standard devia- tion of 1.44 among all the model combinations (Table 4).
However, since K = 0.85, this model combination still overestimated the observed runups and did not meet the requirements for a valid tsunami source based on the Japan Society of Civil Engineers (2002).To address this issue, we performed a simplification for the model combination in this study by reducing the slip amount from the fault model of Imamura et al. (1997).Finally, using a slip amount of 6.4 m for this study combination model, we obtained K and κ values of 1.00 and 1.44, respectively (Table 4).Figure 8 shows the detailed scheme of our study combination model.Moreover, to prove the eyewitness accounts that SL 4 arrived at the southern coast earlier than the main tsunami, and the detailed for the calculated tsunami waveform from our study combination model can be seen in Fig. 9.

Conclusions
We agree with the hypothesis proposed by previous researchers that the 1996 Biak, Indonesia earthquake triggered several SLs around Biak Island.Our analysis revealed two interesting points: 1.The earthquake triggered a large SL 2 and a small SL 4, and the tsunamis generated by these SLs reached the southwest and southern coasts of Biak Island,  respectively, together faster than the main tsunami waves from the earthquake fault.2. During this event, eyewitnesses on the southern coast might have seen the small tsunami wave generated from SL 4, but the observed runups in the coast were due to the main fault earthquake.
Since the models used in previous studies could not explain the runup distribution caused by the SL, we built a model by combining and modifying a previous fault model, SL 2, and SL 4. We chose the fault-based model reported by Imamura et al. (1997), because it yielded a better κ value than the other combined models; however, this combination model did not fulfill the Japan Society of Civil Engineers (2002) criteria for a valid tsunami source.Hence, for this study combination model, we modified the amount of slip of the fault-based model proposed by Imamura et al. (1997) from 8 m to 6.4 m and achieved K = 1.00 and κ = 1.44.
Our study model generally showed good agreement with the observed runups but still cannot reproduce the maximum runup of 7.7 m and the eyewitness tsunami arrival time of 1 min.This limitation may be attributed to several factors: 1.Our multibeam data covered only half of the possible SL area southwest of Biak Island.Therefore, we do not know if another SL candidate was present in the area that was not covered in the data.
2. The observed runup data had inaccurate coordinates.
3. The eyewitness tsunami arrival time had many uncertainties, because the data were obtained from eyewitness accounts, and the coordinates from which they saw the tsunami were unknown.
However, we can attempt to prove the hypothesis even with the data obtained, and the results were more accurate than the results of the previous fault-only models.To maintain consistency and achieve better results, this study should be extended if more complete and accurate observational data become available in the future.

Fig. 1
Fig. 1 Tectonic setting around Biak Island.The black lines indicate local faults.The black lines with triangles indicate subduction zones.The yellow stars indicate the epicenters of tsunamigenic earthquakes.Small colored circles indicate earthquake epicenters.The focal mechanism of the 1996 event is from the Global Centroid Moment Tensor(Ekström et al. 2012)

Fig. 2
Fig. 2 Map of the study area.The yellow star indicates the epicenter of the 1996 Biak tsunamigenic earthquake.Numbers along the coast indicate eyewitness tsunami arrival time in minutes.The black triangles indicate eyewitness observation points by the local people and are used as virtual tidegauges.The red, green, and blue circles indicate the observed runup heights.Layers 1 and 2 correspond to the regions outside and inside the blue rectangle, respectively, on the inset map, and Layers 3a and 3b correspond the interior of the red rectangles

Fig. 3
Fig. 3 Fault models from previous studies with their calculated runup.a Fault models from the previous studies; the yellow star indicates the epicenter of the 1996 tsunamigenic earthquake.b-d Comparisons of the observed and calculated tsunami runups for the northeastern, southwestern, and southern coasts, as shown in Fig. 2, respectively

Fig. 4
Fig. 4 Backward tsunami travel time (TTT) results.The triangles indicate eyewitness observation points, and solid colored lines indicate the backward TTT.The red solid and dashed line circles indicate possible areas and + 1 min buffer possible area of SL, respectively.The numbers indicate TTTs in minutes.The yellow, red, blue, green, and purple rectangles are SL 1, SL 2, SL 3, SL 4, and SL 5, respectively

Fig. 5
Fig. 5 Detailed location and cross-section of the submarine landslides (SLs).a Southwestern SLs, b southern SLs, c detailed cross-section for all SLs.The yellow, red, orange, green, and purple rectangles and ellipses correspond to SL 1, SL 2, SL 3, SL 4, and SL 5 and their SL deposits, respectively.The green solid and dashed circles show the possible area and + 1 buffer of the possible area, respectively, obtained from backward TTT.The black solid line indicates bathymetry depth.The yellow, white, and red dashed lines indicate the Biak Fault Zone, SL traces, and assumed sliding SL planes, respectively

Fig. 6
Fig. 6 Comparison of the calculated tsunami runups in Layers 3-Layer 2. a Layer 3a-Layer 2 and b Layer 3b-Layer 2. The colored circles indicate the residual results of the fault models

Fig. 7
Fig. 7 Comparison of the observed and calculated tsunami runups from SL models only and fault-SL models.a-c Calculated runups from SL models only.d-f Calculated runups from the combination of fault models, SL 2, and SL 4

Fig. 8
Fig. 8 Combination model for the 1996 Biak, Indonesia, tsunami.The purple and red rectangles indicate the modified fault model based on the model proposed by Imamura et al. (1997), SL 2, and SL 4, respectively

Table 2
SL parameters

Table 3
Comparison of tsunami arrival time.Eyewitness tsunami arrival time, calculated tsunami arrival time and first wave impulse at northeastern, southwestern and southern observation points from SL simulation

Table 4 K
and κ results for combination of fault, SL 2, and SL 4 models