Global optimization of a numerical two-layer model using observed data: a case study of the 2018 Sunda Strait tsunami

Following the eruption of Mount Anak Krakatau, a considerable landslide occurred on the southwestern part of the volcano and, upon entering the sea, generated a large tsunami within the Sunda Strait, Indonesia, on December 22, 2018. This tsunami traveled ~ 5 km across the strait basin and inundated the shorelines of Sumatra and Java with a vertical runup reaching 13 m. Following the event, observed field data, GPS measurements of the inundation, and multibeam echo soundings of the bathymetry within the strait were collected and publicly provided. Using this dataset, numerical modeling of the tsunami was conducted using the two-layer (soil and water) TUNAMI-N2 model based on a combination of landslide sources and bathymetry data. The two-layer model was implemented to nest the grid system using the finest grid size of 20 m. To constrain the unknown landslide parameters, the differential evolution (DE) global optimization algorithm was applied, which resulted in a parameter set that minimized the deviation from the measured bathymetry after the event. The DE global optimization procedure was effective at determining the landslide parameters for the model with the minimum deviation from the measured seafloor. The lowest deviation from the measured bathymetry was obtained for the best-fitting parameters: a maximum landslide thickness of 301.2 m and a landslide time of 10.8 min. The landslide volume of 0.182 km3 estimated by the best-fitting parameters shows that the tsunami flow depth could have reached 3–10 m along the shore with a K value of 0.89, although the simulated flow depths were underestimated in comparison with the observation data. According to the waveforms, the general wave pattern was well reproduced at tide gauges during the event. A large number of objective function evaluations were necessary to locate the minimum with the DE procedure to fix the grid cell size to 20 m; this limited the accuracy of the obtained parameter values for the two-layer model. Moreover, considering the generalizations in the modeling of landslide movements, the impact landslide time and thickness must be carefully calculated to obtain a suitable accuracy.


Introduction
Following the eruption of Mount Anak Krakatau, a considerable landslide occurred on the southwestern part of the volcano. This landslide generated a large tsunami around the Sunda Strait, Indonesia, on December 22, 2018, and upon entering the sea, this landslide caused an immense tsunami that traveled ~ 5 km across the strait and inundated the shores of Sumatra and Java with a vertical runup reaching 13 m . Following the event, observed field data, GPS measurements of the inundation, and multibeam echo soundings of the bathymetry within the strait were collected and provided to the public Syamsidik et al. 2020). Figure 1 shows a map of the Sunda Strait in addition to the location of Anak Krakatau in the middle of the strait.
Volcanic eruption-triggered tsunamis have been recorded several times; one of the most recent marine Open Access *Correspondence: pakoksung@irides.tohoku.ac.jp 1 International Research Institute of Disaster Science (IRIDeS), Tohoku University, Sendai 980-0845, Japan Full list of author information is available at the end of the article caldera eruptions was the 1883 Krakatau event. The resulting tsunami caused approximately 35,000 casualties (Giachetti et al. 2012) and affected the coastal area around the Sunda Strait and the western side of the Indian Ocean basin approximately 3000 km from Krakatau. The event produced approximately 20 km 3 of pyroclastic deposits and generated a tsunami runup of 42 m (Choi et al. 2003). The event was modeled by Maeno and Imamura (2011) using the nonlinear shallow water equation, and the 1883 tsunami was successfully reproduced.
The numerical modeling of tsunamis generated by subaerial-submarine landslides requires different methods if the conclusions are displayed as completely coupled schemes. Numerical models comprising coupled dynamic schemes for landslide-generated tsunamis are limited (Heinrich et al. 2001;Shigihara et al. 2006), while those for earthquake-generated tsunamis are developed and widely implemented. There have been several studies of landslide-generated tsunamis using nonlinear shallow-water hydrostatic models (Kowalik and Murty 1993;Satake 1995;Heidarzadeh et al. 2014), Boussinesq nonhydrostatic models ), the MOST model (Titov and González 1997), the Cornell Multi-grid Coupled Tsunami Model (COMCOT) (Liu et al. 1998), and Tohoku University's Numerical Analysis Model for Investigation of Nearfield Tsunamis, No. 2 (TUNAMI-N2) model (Imamura and Imteaz 1995). Moreover, several studies have been performed to model the tsunami-generation mechanisms using other uncoupled methods, in which the generation and propagation periods are episodic (Lastras et al. 2005;Iglesias et al. 2012). One such example is the conceptual model of the BIG'95 submarine landslide and the numerical simulation of the propagation of the generated tsunami (Skvortsov 2002;Watts et al. 2003;Fine et al. 2005;Tinti et al. 2006Tinti et al. , 2011. In the case of the Palu tsunami, the submarine landslide source was implemented by using a combination of models, namely, Titan2D (Pitman et al. 2003;Patra et al. 2005;Titan2D 2016) and JAGURS (Baba et al. 2017), with a multi-landslide location in the bay. The volume of the multi-landslide source varied from 0.02 to 0.07 km 3 among six different locations (Nakata et al. 2020). In addition, the TUNAMI-N2 model was applied to model the submarine landslide-induced tsunami in Palu Bay by Pakoksung et al. (2019), who proposed that the main source of the landslide that generated the 2018 Palu tsunami was in the northern part of the bay. This Anak Krakatau tsunami event was recently modeled (Grilli et al. 2019) using the 3D Nonhydrostatic Wave (NHWAVE) model (Ma et al. 2012;Kirby et al. 2016) to simulate the landslide and the 2D Boussinesq FUNWAVE-TVD model  to simulate the tsunami propagation. Furthermore, COMCOT (Liu et al. 1998;Wang and Liu 2006) was also used to model this event (Heidarzadeh et al. 2020), and the initial sea surface elevation was identified as the landslide source. Some previous studies implemented other methods to propose the landslide source for the 2018 Anak Krakatau tsunami. For example, synthetic aperture radar and broadband seismic observations were used to identify the landslide volume (< ~ 0.2 km 3 ), as presented by Ye et al. (2020). In this study, the two-layer TUNAMI-N2 model (Imamura and Imteaz 1995;Pakoksung et al. 2019) is implemented to model the 2018 Anak Krakatau tsunami.
The 2018 Anak Krakatau tsunami provides a specific opportunity to investigate and capture the essential mechanism responsible for the generation and propagation of a landslide-induced tsunami through numerical modeling, as an exhaustive amount of data was acquired to calibrate the model. Optimization was achieved by comparing the field observations and deposit simulation results to delineate the impact physics of tsunamis without the scale effects detected by the physical experiment.
Models are usually calibrated by generating a set of parameters to match real data. Having a good model and a strong reliable numerical method for solving problems is as important as performing a good parameter adjustment of the model according to physical measures. On the other hand, a good model together with a good numerical method can lead to totally incorrect results with poorly calibrated parameters. Only a few studies have applied the global optimization method to tsunami modeling for rockslide-induced tsunamis, as proposed by Gylfadottir et al. (2017).
Our work is similar to the recent study by Gylfadottir et al. (2017) in that we applied a global minimization algorithm to calibrate a rockslide-induced tsunami model. The approach of Gylfadottir et al. (2017) consists of minimizing, with a differential evolution (DE) procedure, a cost function corresponding to the modelobservation good-fit of a rockslide model to generate a tsunami in Lake Askja. The similarities between this work and our research are twofold. First, a boundconstrained global stochastic minimization algorithm is used. Second, the method to assess the optimality of the achieved solution uses a pool of independent and randomly initialized minimization experiments. However, the method we are recommending contrasts from their mechanism. Based on the different natures of the rockslide model and the two-layer model, the two-layer tsunami model is based on the flow of two different fluids, while the rockslide model is represented by solid block movement.
The aim of this study is to investigate the possible source of the 2018 Sunda Strait tsunami using preliminary data. This study focuses on a tsunami generated by a subaerial-submarine landslide. The method employed to find the tsunami source is to apply a global optimization procedure to determine the unknown landslide parameters through a comparison between the measured and simulated bathymetry with an inverse modeling methodology, using a two-layer model.
The remainder of this paper is presented as follows. The methodology is explained in "Methodology" section. The beginning of "Methodology" section presents the flank collapse of the Anak Krakatau volcano for the studied event ("Anak Krakatau flank collapse event in 2018" section). Then, the topography/bathymetry change from the event is described in "Topography and bathymetry change from the 2018 Anak Krakatau flank collapse event" section, and the modeling framework of global optimization combined with the two-layer model is explained in "Modeling framework" section. "Global optimization tool" section presents the optimization algorithm, and the two-layer model and the setup of the model parameters are described in "Twolayer model" section. Finally, the modeling results are presented in "Results" section, the results are discussed in "Discussion" section, and the conclusions are provided in "Conclusions" section.

Anak Krakatau flank collapse event in 2018
The flank collapse of the Anak Krakatau volcano in 2018 generated a subaerial-submarine landslide tsunami in the Sunda Strait as a result of the movement of soil materials originally located on the upper slope. This process produced a subaerial-submarine landslide on the land and seafloor, forming a submarine deposit. The tremendous landslide was probably triggered by the volcanic eruption of Anak Krakatau. This subaerial-submarine landslide should have produced a tsunami that affected the coastal area of the Sunda Strait. The shoreline change on Anak Krakatau Volcano Island is shown in Fig. 2a (before) and Fig. 2b (after); these changes can be detected by satellite images from before and after the eruption. These flank collapse landslides, which were located in different areas of the coastline, were the main cause of the studied tsunami.

Topography and bathymetry change from the 2018 Anak Krakatau flank collapse event
Complete bathymetric and topographic data from throughout the Sunda Strait and surrounding continental areas were provided by BATNAS and DEMNAS, Indonesia (available at https ://tides .big.go.id/DEMNA S/ index .html), respectively. The data consist of two datasets with resolutions of 180 m for the bathymetry in the sea and 8 m for the topography on land and represent the bathymetry and topography before the studied event. Both datasets were resampled to three domains with resolutions of 20, 60, and 180 m, as shown by the location of 20 m in Fig. 1, while the 60 m area was buffered by approximately 1 km from the border of 20 m. The domain of 180 m was directly used, and the data were only cropped to cover the Sunda Strait area. The 20-m and 60-m resolutions were generated by using the Kriging method in QGIS, based on the 8 m data for the land and the 180 m data for the sea. Only the 20 m (orange box in Fig. 1) around the Anak Krakatau area was used to run the two-layer model, while the other 20 m (blue box in Fig. 1) was used to simulate the tsunami runup in the observed area.
The Naval Hydrographic and Oceanographic Center of Indonesia used multibeam sonar equipment to survey the bathymetry after the eruption of Anak Krakatau (PUSH-IDROSAL-Naval Hydrographic and Oceanographic Center 2019) as proposed by Muhari et al. (2019). The measured area of new bathymetry data covers only the front of the collapse caldera, as shown in Fig. 3, and the observed area covers only the southwestern part of Anak Krakatau with an average difference of 30 m from the bathymetry before the eruption. Based on the difference between both bathymetry datasets, we question when the collapse of soil materials from the volcano eruption reproduced the shape approximating the observation data (the bathymetry after the event).

Modeling framework
The general framework of this study, which is shown in Fig. 4, indicates that the main input data required by the two-layer model are topography/bathymetry data and landslide data (shape and depth). Topographic/bathymetry data were resampled from the DEMNAS and BAT-NAS databases, Indonesia, which provide a dataset of geographic information. The landslide geometry was obtained from the satellite images before and after the event. The two-layer model was integrated with a global optimization tool (DE) to calibrate the sensitive parameters of the landslide. The output results of this model were landslide deposition and tsunami propagation. The landslide deposition results were verified by using the bathymetry at the area around the flank collapse observed using multibeam sonar equipment.

Global optimization tool
Based on the uncertainties in the parameters of landslide movement in the two-layer model, an inverse method was applied to determine the optimal parameters that best describe the bathymetry after the volcanic eruption. For the two-layer model, the best-fitting parameters were estimated by a grid search, whereby different combinations of parameters were repeatedly run on the model. A set of initial models was used to obtain an approximate idea of the landslide parameter ranges. The objective function, namely the quality of fit, was considered by using the sum squared error (SE) index as a summation of the squared differences between the observed and simulated bathymetry values in front of the volcano.
The grid search method initially considers the result of the SE index in a parameter to be complex with a gradient-based minimization process, and the algorithm was regarded as unsuccessful if the solution became stuck in a local minimum. Then, a global optimization procedure was selected to test the parameter space. The DE procedure provided by Storn and Price (1997) was selected and applied to the two-layer model, as shown by the flowchart in Fig. 5. The DE approach is a simple method and uses a few control variables to control the minimization. This method has good performance compared to other minimization procedures, such as the genetic algorithm (GA) and particle swarm optimization (PSO) methods. DE is similar to the GA in the sense that it uses the same evolutionary operators such as mutation, crossover, and

Two-layer model
The tsunami model in the Sunda Strait was assessed with a two-layer numerical model that was developed to solve the nonlinear shallow water equation with two interfacing layers and appropriate kinematic and dynamic boundary conditions at the land and seafloor, their interface, and the water surface (Imamura and Imteaz 1995;Pakoksung et al. 2019). This two-layer model simulates landslide-generated tsunamis by modeling the interactions between the generated tsunami and subaerialsubmarine landslides as the upper and lower layers, respectively. Subaerial-submarine landslides produce tsunamis similar to those produced by earthquakes with a vertical displacement of the seafloor creating a similar displacement at the sea surface, as shown in Fig. 6. The important point of the model is the landslide thickness, which is used to identify the boundary of landslide movement.
The mathematical model employed in the landslide tsunami code consists of a stratified medium with two layers, as shown in Fig. 6. The first layer is composed of a homogeneous inviscid fluid with a constant density ρ 1 representing seawater, and the second layer is based on a fluidized granular material with a density ρ g and a porosity ϕ . In this study, the mean density of the fluidized debris is constant and equals ρ 2 = (1 − ϕ)ρ g + ϕρ 1 , as noted in previous research (Macías et al. 2015). The two fluids, water and fluidized debris, are hypothesized to be immiscible in this study. The governing equations are written as follows: Continuity equation of the first layer: Momentum equations of the first layer in the X and Y directions: Two-layer conceptual schematic of the parameters for modeling the subaerial-submarine landslide Pakoksung et al. Geosci. Lett. (2020) 7:15 Continuity equation of the second layer: Momentum equations of the second layer in the X and Y directions: where index 1 relates to the upper layer and index 2 indicates the second layer. Z i (x, y, t) , i = 1, 2 is the level of the layer at each point (x, y) at time t , where the level value is measured from a given reference level. Q i (x, y, t) , i = 1, 2 is the vertically integrated discharge in the x and y directions. g is the gravitational acceleration. ρ 1 and ρ 2 are the densities of the first and second layers, respectively. τ i (x, y, t) is the bottom stress in each layer at each point (x, y) at time t . In the momentum equation, the interaction between the first layer and the second layer is determined by the fifth term of the momentum equation. The landslide movement is stopped in the first layer based on the sliding time, which is one of the limitations on the two-layer model used in this study. The parameters were kept fixed as follows: water density = 1000 kg/ m 3 , landslide density = 2,000 kg/m 3 , and Manning coefficient = 0.025.
The hypothesized landslide modeling for the two-layer model of the 2018 Anak Krakatau tsunami is shown in Fig. 2c; these landslides can be detected by satellite images from before and after the eruption. The change in the coastline was used to consider the landslide model using ellipsoid modeling. The shape of the landslide hypothesized by ellipsoid modeling has a length of approximately 1.7 km and a width of approximately 0.8 km; the center is located at a latitude of − 6.102 degrees and a longitude of 105.4°.
The remaining parameters required as input to the two-layer model area sliding time and landslide depth could not be captured by the satellite image and field survey. Hence, optimization modeling was implemented to achieve the optimum values of both parameters. These parameters can be roughly estimated to guide the optimization. The landside depth proposed in this study was based on the results of Giachetti et al. (2012), who mentioned the worst case Anak Krakatau flank collapse with a volume of approximately 0.28 km 3 . To fit the landslide into the landslide area in this study, the landslide depth would have been approximately 450 m. The approximated landslide depth from the worst-case scenario was identified as the upper boundary value for optimization. The lower boundary of the landslide depth was assumed to be one-third of the maximum height of the volcano.
The impact sliding time is difficult to estimate due to the inherent uncertainty in the speed as a result of the frictional parameters. Then, the sliding time of the landslide in this study was assumed to range from 1 to 20 min.

Results
The numerical tsunami simulation results from the possible subaerial-submarine landslide sources were estimated by the global optimization tool (DE method). The original positions and final positions of the sources were investigated to assess the impact of a tsunami generated by each volcanic eruption-induced landslide. The height of the potential tsunami on the coastal area was investigated based on the maximum tsunami amplitude at the shoreline and coastal area defined as the maximum tsunami height from the possible submarine landslide scenario. The maximum tsunami amplitude and inundation area from the numerical modeling were calculated from the terrain data with a grid size of 20 m over a simulation time of 90 min.

Global optimization procedure result
The global optimization method using the introduced DE algorithm was applied to determine the optimal landslide parameters though the two-layer model. that were collected based on the minimum SE of the last iteration. Overall, the fit was very good, and the improvement for the underwater area extended far from the shore (approximately 1 km).

Tsunami simulation resulting from the best-fitting parameters of the landslide
Numerical tsunami simulations were performed based on the best-fit parameters of landslide sources to simulate possible tsunamis and to collect numerical results along the coastline of the Sunda Strait. Three bathymetry and topography domains with resolutions of 20, 60, and 180 m were used to perform a constant-grid tsunami simulation, as shown as the blue box in each domain in Fig. 1. The simulation yielded some 11 million equations and unknowns to be solved at each time step of 0.01 s. At the boundary lines, the open sea was limited with nonreflective boundary conditions, whereas the coastal areas had no specific boundary conditions for wet/dry fronts (Imamura 1995).
The mathematical model parameters were modeled to match the simulated submarine deposits, as shown in Fig. 3, which presents a subaerial-submarine landslide that slid southwestward down a slope of approximately 20°. This landslide movement was modeled from the starting position to the end position, and the sliding is oriented downslope along the sliding plane. Therefore, if the hypothesized submarine landslide occurred, we expect that it would have reproduced the associated tsunami and its impacts on this study area. Simulation snapshots of the propagating tsunami at different times are shown in Fig. 8 for 0, 5, 10, 20, 30, and 60 s after the tsunami was generated by the volcano eruption-induced landslide. The modeled landslide mostly collapsed within approximately 1 min of the volcano altitude decreasing from approximately 300 m to approximately 100 m. The flank collapsed via the slow sliding of a 0.182-km 3 mass down the subaerial-submarine slope, and sliding completely stopped after 10.8 min (as determined by the optimization procedure).
The landslide exerted a barrier effect that is reflected in the maximum amplitude distribution, with the subaerialsubmarine zone in which sliding initiated exhibiting the maximum amplitudes. A region of highly elevated amplitudes can also be observed surrounding the volcanic mountain; however, the amplitudes decrease dramatically in the middle between the volcano and the shorelines of the Sunda Strait in Sumatra and Java. On the other hand, the amplitude increases in the surf zone near the shorelines of both islands because of the shallow water depths therein. The maximum tsunami height from the numerical simulation based on the best-fitting parameters is shown in Fig. 9. The flank collapse generated the maximum positive elevation (approximately 70 m on the southwestern side of the volcano approximately 1 km from the sliding plane).
Furthermore, a velocity decrease was noted over the main high relief, and this velocity variation produced an alteration in the concentric pattern of arrival times obtained from the best-fitting parameters covering the domain area, as shown in Fig. 10. The distribution of arrival times appears marked by inflection lines in the propagation pattern. In this regard, the presence of continuous seamounts in the middle of the strait amplified the amplitude; however, these seamounts were responsible for wave deceleration and therefore produced a delayed arrival along the coastline. The first place struck by the tsunami was Surtung Island (far from the volcano, approximately 3 km). Moreover, Rakata Island was affected at the same time (approximately 60 s). On Java, the first tsunami arrived after approximately 35 min at Pandeglang city (approximately 51 km from the volcano). The first tsunami arrived at Sumatra Island in the city of  The effects of the tsunami (from the best-fitting parameters) along the coastline of the Sunda Strait are considered in the three areas mentioned in Fig. 1, as shown in Figs. 12, 13, and 14, based on the surveyed area in Syamsidik et al. (2020). In the flood area along the coastal zone in area A, the tsunami flow depth has an average value of approximately 2.75 m with a maximum of 5.8 m in the middle of this area because of the large island (Sangiang Island) situated at the entrance of the strait. A comparison between the simulated wave height with the surveyed result revealed an underestimation, as shown in Fig. 12. Figure 12a presents the distribution of survey points and the computed flood areas along the coastline, indicating that all of the surveyed points were located in the computed flood area. Profiles of the observed and simulated flow depths along the coastline are compared in Fig. 12b, which reveals that the simulated results are close to the observed data in the northern part of the strait, while underestimated flow depths are observed in the south. In area B, the tsunami flow depth has an average value of approximately 3.0 m with a maximum of 6.2 m in the north. This area was highly impacted by the tsunami because of the large island (Panaitan Island) in the south, which reflected the tsunami to flow toward this area. A comparison between the simulated and observed flooding in this area (shown in Fig. 13) suggests that the model can simulate a flood covering all the survey points presented in Fig. 13a. Profiles of the observed and simulated flow depths along the coastline reveal that the simulated results are similar to the observed data in the north and south, while the flow depths are underestimated is in the middle, as shown in Fig. 13b. In area C, the tsunami flow depth has an average value of approximately 3.3 m with a maximum of 4.5 m in the east. The area is close to the volcano (approximately 80 km) but contains two large islands (Sebuku Island and Sebesi Island) located in the front as natural protection; consequently, this area was impacted less than areas A and B. Figure 14a illustrates that the model can simulate a flood covering all the survey points. Profiles of the observed and simulated flow depths along the coastline reveal that the simulated results are close to the observed data in the west, while the flow depths are overestimated in the east, as shown in Fig. 14b.   Fig. 9 Maximum wave height in the whole Sunda Strait area based on the optimum parameters obtained from the optimization modeling Pakoksung et al. Geosci. Lett. (2020) 7:15

Discussion
The tsunami water levels and flow depths were collected from a field survey after the 2018 Anak Krakatau tsunami by Syamsidik et al. (2020). The tsunami runup was further determined using a numerical model since insufficient field measurement data were recorded. The numerical tsunami model was verified using a performance parameter, such as the geometric mean K or geometric standard deviation κ (Aida, 1978). The K value, which refers to the deviation or variance from the proportion between observed and simulated data, was derived from the mean of the K and κ values. The Japan Society of Civil Engineering (JSCE) recommends values of 0.95 < K < 1.05 and κ < 1.45 for a model to achieve good agreement when evaluating the source of a tsunami and its propagation model (Japan Society of Civil Engineering (2002); Suppasri et al. 2011;. The parameters K and κ are determined as shown below: where x i and y i are the observed and simulated data, respectively, at point i.
The tsunami flow depths from the numerical simulation based on the best-fitting landslide model (achieved by the optimization procedure) and the bathymetric data were evaluated using field measurement data. The tsunami flow depths were compared with the measurements

Fig. 10
Arrival times of the leading tsunami measured by the positive amplitude Pakoksung et al. Geosci. Lett. (2020) 7:15 provided in a previous study (Syamsidik et al. 2020). The comparison reveals that K is approximately 0.89 and κ is approximately 1.70, which shows fair agreement with the above-mentioned JSCE standard due to the effect of building occupation on tsunami runup modeling, as explained in the next paragraph. The observed tsunami flow depth varied from 0.5 to 6.5 m, with the highest values located in the southern part of the beach in Pantai Cisiih. Scatter plots of the observed and simulated data are shown in Fig. 15, based on which the simulation results contain an error relative to the model in the range of approximately 3 m as the standard error. This error is The K and κ values do not agree with the JSCE standard because the governing equation of tsunami runup ignores the effect of a building inside the model, although the building obstructs the flow on the terrain. The obstruction of a building can increase the flow resistance relative to a model without the building, and therefore, the simulation of flow through a building can raise the water level close to the real situation (Copeland 2000;Dutta et al. 2007;Aburaya and Imamura 2002;Fukui et al. 2019). The inclusion of buildings within the tsunami runup modeling process constitutes the limitation of this study; hence, we suggest inputting buildings into the model in future work, which might provide flow depths close to the observation data.
This study presents a preliminary model of the flank collapse associated with the 2018 Anak Krakatau landslide, which generated a tsunami, and the landslide parameters are optimized based on the observed bathymetry after the event. The main purpose of this study is to thoroughly understand the landslide that occurred and estimate its volume despite lacking the submarine landslide mechanism. This study proposes a collapse volume of 0.182 km 3 to obtain the flow depth distribution along the coastline around the Sunda Strait with a K of approximately 0.89 for a comparison with observation data. For comparison with previous studies, our landslide volume  Pakoksung et al. Geosci. Lett. (2020) 7:15 is different from that (approximately 0.27 km 3 ) reported by Grilli et al. (2019). On the other hand, our landslide volume is close to that of Paris et al. (2020) with 0.15 km 3 . Furthermore, our landslide volume is in the range (< ~ 0.2 km 3 ) suggested by Ye et al. (2020), who estimated the volume by using different methods, that is, a combination of synthetic aperture radar and broadband seismic observations. Regarding the four tide gauges, our results are similar to the findings of Grilli et al. (2019) and Paris et al. (2020) considering the pattern of the waveform. Additionally, our maximum wave amplitude (approximately 70 m) is close to the maximum of approximately 80 m proposed by Paris et al. (2020); in contrast, our maximum wave height is quite different from that (approximately 100-150 m) proposed by Heidarzadeh et al. (2020). The best-fitting parameters, namely, a maximum landslide thickness of 301.2 m and a landslide time of 10.8 min with a generated landslide volume of 0.182 km 3 , can simulate the waveform pattern well; accordingly, the general wave pattern was well reproduced at the tide gauges during the event. Our estimated landslide volume is different from that of Grilli et al. (2019) but close to that of Paris et al. (2020) because of the model dimensions. Moreover, the landslide modeling conducted by Grilli et al. (2019) was achieved by using a 3D non-hydrostatic model, while Paris et al. (2020) used 2D modeling to simulate the landslide. A 2D model is also implemented in this study, but our model has a different soil property term (friction angle of sliding material), which was excluded from the governing equation employed in this study. This missing term is a limitation of this study insomuch that the modeling is accomplished without the friction angle of sliding. This constitutes a small difference in the maximum amplitude of the wave (approximately 10 m) from Paris et al. (2020) as a result of the Boussinesq modeling and the resolution of the bathymetry data.
Even though the results are sufficiently good to facilitate a comparison with observation data, some limitations of this study should be explained. First, the landslide movement was produced by several processes, but this study considers only three parameters: the landslide density is fixed, while the other two parameters (depth and sliding time) are varied. These two parameters were selected through the optimization process by minimizing a different deviation from the bathymetry after the event; this selection might have affected the accuracy of the model. In addition, in future work, we suggest that some soil parameters, such as the friction angle, which depends on the soil type, and the friction angle and roughness coefficient between the landslide and sliding plane, be added to the two-layer model (system 1 in this study). Adding soil parameters into the two-layer model was already recommended in previous studies (Ioki et al. 2019;Iverson and Denlinger 2001).  Fig. 15 Scatter plot between the observed and simulated data (from the best-fitting parameters) for the maximum flow depth Pakoksung et al. Geosci. Lett. (2020) 7:15 Conclusions A large landslide released following the southwestern flank collapse of the Anak Krakatau volcano was modeled by a two-layer model integrated with a global optimization model to simulate the subsequent tsunami in the Sunda Strait. The parameters of the landslide model were estimated by global optimization using the bathymetry surveyed after the event. The lowest deviation from the measured bathymetry was obtained to determine the bestfitting parameters: a maximum landslide depth of 301.2 m and a landslide time of 10.8 min with a generated landslide volume of 0.182 km 3 . The results show that the tsunami flow depth reached up to 1-13 m along the shore with a K of 0.89 with only minor underestimation compared with the observation data (flow depth). According to the waveforms, the general wave pattern was well reproduced at the tide gauges during the event. However, the model is restricted by some limitations. The landslide model was hypothesized by considering satellite images for the shape and employing ellipsoid modeling for the volume. Furthermore, the twolayer model was constructed without soil parameters, such as the friction angle and friction roughness. However, the tsunami simulation results related to this landslide are quite consistent with observed flow depths.
The two-layer model proposed in this paper produces a realistic situation that agrees with the main characteristics of the available tsunami data. However, future work should improve our model as new geological data from the volcano becomes available. Compared with the observation data, the results from the two-layer model, such as the flow depth and bathymetric change, were underestimated. The model requires high-resolution terrain data to produce a high-accuracy result and requires the integration of a non-hydrostatic term into the governing equation (Maeda et al. 2016;Gylfadottir et al. 2017). Hence, future efforts should be directed to producing more realistic estimates by modifying the model. Using the sliding time to stop the landslide is the main limitation of this study because the two-layer model does not include a soil property in the governing equation (Eqs. 3,4). Normally, landslides are modeled by adding soil properties such as the Coulomb friction angle or Basal friction angle into the momentum equation of the landslide model (see Harbitz et al. 2003;Paris et al. 2020).
Although there are many uncertainties in tsunami hazard evaluation, such as uncertainties in the submarine landslide geometry, bathymetric and topographic datasets, and numerical modeling, as mentioned above, our intention was not to perform a probabilistic hazard evaluation. Instead, this study proposes that a set of subaerial-submarine landslides triggered by the 2018 volcanic eruption of Anak Krakatau could have produced the Sunda Strait tsunami.