Modelling 2018 Anak Krakatoa Flank Collapse and Tsunami: Effect of Landslide Failure Mechanism and Dynamics on Tsunami Generation

The 2018 Anak Krakatoa volcano flank collapse generated a tsunami that impacted the Sunda Strait coastlines. In the absence of a tsunami early warning system, it caused several hundred fatalities. There are ongoing discussions to understand how the failure mechanism of this event affected landslide dynamics and tsunami generation. In this paper, the sensitivity to different failure scenarios on the tsunami generation is investigated through numerical modelling. To this end, the rate of mass release, the landslide volume, the material yield strength, and orientation of the landslide failure plane are varied to shed light on the failure mechanism, landslide evolution, and tsunami generation. We model the landslide dynamics using the depth-averaged viscoplastic flow model BingClaw, coupled with depth-averaged long wave and shallow water type models to simulated tsunami propagation. We are able to match fairly well the observed tsunami surface elevation amplitudes and inundation heights in selected area with the numerical simulations. However, as observed by other authors, discrepancies in simulated and observed arrival times for some of the offshore gauges are found, which raises questions related to the accuracy of the available bathymetry. For this purpose, further sensitivity studies changing the bathymetric depth near the source area are carried out. With this alteration we are also able to match better the arrival times of the waves.


Introduction
Landslides and volcano flank collapse tsunamis arguably represent the second most important tsunami source with respect the hazard they pose to populations Tappin 2010). From large sub-aqueous landslides, involving volumes of hundred to thousands km 3 (Masson et al. 2006), to smaller subaerial landslides impacting the water surface with large speed , landslides span very different scales and types of generation mechanisms (Heidarzadeh et al. 2014;Løvholt et al. 2015; Yavari-Ramshe and Ataie-Ashtiani 2016). Volcano flank collapses comprise one of the most important subset of tsunamigenic landslides, as recently demonstrated by the 2018 Anak Krakatoa tsunami (Babu and Kumar 2019;Gouhier and Paris 2019;Grilli et al. 2019;Walter et al. 2019;Muhari et al. 2019;Paris et al. 2020;Heidarzadeh et al. 2020;Novellino et al. 2020;Ye et al. 2020;Putra et al. 2020) and is the focus of this study. History has shown that Anak Krakatoa, by no means, represents a unique event.
In the seventeenth century, two devastating flank collapse tsunamis, in Shimabara 1792 (Sassa et al. 2016), and in Oshima-Oshima 1741 (Satake and Kato 2001;Satake 2007) impacted Japanese coastlines. Both caused several thousand fatalities, the 1792 Shimabara event being the most fatal of all historical landslide tsunamis. It has also been speculated that flank collapses were a major factor in generating the 1883 Krakatoa tsunami (Nomanbhoy and Satake 1995). Reviews of historical events in Indonesia ) and the Caribbean  list several more instances of landslide tsunamis from volcanoes. In 2002, a smaller event emerged from the Stromboli volcano (Tinti et al. 2005). Consequently, a more widespread tsunami threat emerging from potentially much more voluminous flank collapses with oceanic propagation have been subject to speculation (Ward and Day 2001). To test this hypothesis, this had lead to development of models for dealing with both the near-field and far-field tsunami generation and propagation Abadie et al. 2012). Subject to debate has been to which extent flank large collapses occur as single events, or with stage-wise release of the volume as evident from deposits of Canary Island events (Wynn et al. 2002;Masson et al. 2006;Paris et al. 2017). With the recent Anak Krakatoa event in mind, there is a timely need to better understand how this mechanism controls the tsunami generation.

The 2018 Anak Krakatoa Event
On 22 December 2018, the collapse of Anak Krakatoa produced a devasting tsunami that inundated southern Sumatra and western Java in Indonesia, and caused several hundred fatalities (Muhari et al. 2019;Putra et al. 2020). The tsunami hit the nearest coast 35 min after the southwestern volcanic flank collapse of Anak Krakatoa. Four tidal gauges recorded the tsunami heights at the western Java coast, in Marina Jambu and Ciwandan, and at the southern Sumatra coast, in Panjang and Kota Agung (see Fig. 1). The maximum surface elevation of these four gauge stations was recorded in Marina Jambu with a height of 1:40 m. A post-tsunami field survey by Muhari et al. (2019) in two regions at the western Java coast found a maximum measured runup height above 13 m in Pantai Legon Tanjungjaya, south of the western Java coast (see Fig. 1). Putra et al. (2020) reported a similar maximum runup height in the same region.
The Krakatoa complex, containing four islands, the Anak Krakatoa Volcano, Sertung Island, Panjang Island, and Rakata Island, is located in the centre of Sunda Strait, Indonesia. The three surrounding islands, together with the submarine caldera, are remnants of the 1883 Krakatoa eruption, which has been studied by several authors (Francis and Self 1983;Camus and Vincent 1983;Decker and Hadikusumo 1961;Sigurdsson et al. 1991;Deplus et al. 1995;Nomanbhoy and Satake 1995;. Anak Krakatoa has frequently erupted since 1927, with eruptions typically strombolian to vulcanian in style, characterized by small explosive eruptions with columns reaching to 1 km in height, pyroclastic and lava flows (Camus et al. 1987;Kimata et al. 2012;Suwarsono et al. 2019). During these eruptions in the last decades, Anak Krakatoa built up to a height of around 330 m above sea level. Based on seismographs on the Sertung island and elsewhere on Sumatra and Java, Muhari et al. (2019) identified the 2018 Anak Krakatoa flank collapse to have taken place at 20:55 local time. Giachetti et al. (2012) simulated a hypothetical Anak Krakatoa flank collapse already several years ago, and their scenario volume (0:28 km 3 ) and configuration was in fact very similar to the event that occurred in 2018, but their simulated waves were substantially larger. Grilli et al. (2019) modelled the 2018 Anak Krakatoa flank collapse as a deformable landslide on a geometrical failure plane interpreted from their satellite images. They applied three different volumes and two rheologies, one being a granular material and the other one being a dense viscous fluid. Moreover, Paris et al. (2020) used a Savage Hutter model to simulate the landslide dynamics and coupling to the tsunami generation. A sensitivity study on the landslide friction angle was carried out, and a reasonably good agreement with surface elevations observations at the four gauge stations were found. A moderate effect of frequency dispersion for the wave propagation was also reported. Heidarzadeh et al. (2020) studied a combination of numerical simulations with physical sliding experiments to propose an initial wave of the event. The influence of the geometry of the initial wave and related energy were tested. However, the degree of sensitivity studies in above studies were limited to study the variability to landslide volume and model rheologies (Grilli et al. 2019), to the friction angle (Paris et al. 2020), to the geometry of the initial wave (Heidarzadeh et al. 2020), and none of the above studies investigated the sensitivity to the rate of landslide mass release, bathymetric depth, failure orientation, or inundation, which are all subject to investigation here.
In this study, we treat the 2018 Anak Krakatoa flank collapse as a depth-averaged visco-plastic Herschel-Bulkley type flow with the landslide model BingClaw (Løvholt et al. 2017;Kim et al. 2019;Vanneste et al. 2019). BingClaw has previously been successfully used for modelling clay-rich materials of the largest submarine landslides of the world such as Traenadjupet and Storegga (Løvholt et al. 2017;Kim et al. 2019), and the 1929 Grand Banks landslide (Løvholt et al. 2018;Zengaffinen et al. 2020). Here we use BingClaw for modelling the volcanic rockslope landslide failure dynamics. The reason why we choose to use BingClaw for this study is that it is the only model we are aware of that allows the failure mass to be released gradually over time, rather than as a single failure.

Aims and Overview
The overarching aim of our study is to understand the mechanism of the landslide that caused the tsunami. In order to do so, we address the following questions.
1. Was the landslide a single stage or a multistage event? 2. How does the landslide directivity affect tsunami generation? 3. How sensitive is the model to potential errors in the bathymetry?
Section 2 presents the landslide and tsunami models, gives an overview of relevant input parameters and modelling setups, and describes the tidal gauge and runup height observations. In Sect. 3 we use a basic scenario to describe model tests, and to show contour plots of the landslide and tsunami. Then, through a parametric sensitivity analysis, we show the effect of the model parameters on the surface elevations at the four main gauge stations and compare the results to observations. For the basic scenario, which is the best fit scenario, an inundation study is compared to field observations. Section 4 summarises key conclusions.

Landslide Model
In this paper, we use the viscoplastic landslide model BingClaw (Løvholt et al. 2017;Kim et al. 2019;Vanneste et al. 2019) to simulate the volcano flank collapse of the 2018 Anak Krakatoa event. The model implements the Herschel-Bulkley rheology in a two-layer depth-averaged formulation (Huang and Garcia 1998). Kim et al. (2019) presented a detailed description and derivation of the model. BingClaw solves the mass conservation equation integrated over the entire flow depth, the momentum conservation equation integrated separately over the plug layer depth and shear layer depth, in two horizontal dimensions (2HD). The model combines a finite volume method for the leading order terms with a finite difference method for the source terms, and employs the conservation law package ClawPack (Mandli et al. 2016) using the GeoClaw module ). If the earth pressure p ¼ q d g 0 hrðh þ bÞ exceeds the shear strength of the material in a given computational cell, the dynamic equations are solved with a Godunov fractional step method. First, the equations without friction terms are solved, then the frictional terms. If the shear strength is larger than the earth pressure, no motion is imposed. A notable aspect of the model is that it can simulate remoulding of landslide material, changes in material properties over time. The maximum yield strength of the material will be lowered to a residual yield strength as a function of accumulated shear deformation, which is described by the use of a heuristic formulation (De Blasio et al. 2005) s y ðcÞ ¼ s y r þ ðs y 0 À s y r Þe ÀCc ð1Þ where s y r is the residual yield strength, s y 0 the maximum yield strength, C the remoulding rate, and c the accumulated shear deformation.
In addition to remoulding of the yield strength over the entire landslide area, BingClaw allows the user to pre-remould parts of the landslide volume before the landslide motion is initiated. Here, we preremould a subset of the landslide volume extending some distance up from the slide toe. This feature allows the down-slope part to flow out first, which can help setting up landslide simulations that mimic a cascading failure.

Tsunami Models
We use two depth-averaged long wave models to simulate propagating tsunamis over varying bathymetry. The first model is GeoClaw that is used for the majority of the analyses. It assumes hydrostatic pressure, captures propagation of breaking waves, bottom drag, and dry-land inundation with a moving shoreline Berger et al. 2011). GeoClaw incorporates the non-linear shallow water (NLSW) equations in conservative form and in 2HD, and uses the shock-capturing Finite Volume methods and Riemann solvers for ClawPack. Secondly, we use the GloBouss model, that is run with both linear shallow water (LSW) and linear dispersive (disp) equations in order to evaluate the importance of dispersion. GloBouss lacks features such as a moving shoreline. Hence, it cannot be used for simulating dry-land inundation (Pedersen 2008;Løvholt et al. 2008).
The temporal volumetric change of the bathymetry caused by the submarine mass movements are the primary sources for the tsunami generation. Such progressing changes in the bathymetry are output from BingClaw at given times and are used directly in the tsunami generation model. We use only GeoClaw as a tsunami generation model. GloBouss is merely used as a tsunami propagation model that takes the surface elevations and velocities from GeoClaw at 600 s as initial conditions. At this time the landslide is nearly at rest and the leading waves have escaped the Krakatoa complex.

Geometrical and Geotechnical Input
We employ the default parameters listed in Table 1. Our landslide density q d ¼ 1500 kg m À3 is based on the study of a potential Anak Krakatoa flank failure by Giachetti et al. (2012). The Herschel-Bulkely flow exponent n is 0\n\1 for shear-thinning materials. During preliminary tests, it was found that n is not important for the tsunami generation (Zengaffinen et al. 2020). We employ n ¼ 0:5. The reference strain rate _ c r relates dynamic viscosity, yield strength and the Herschel-Bulkley flow exponent . We set _ c r ¼ 1:6 Á 10 7 s À1 . The maximum yield strength s y 0 of the landslide material defines its strength at the initial time unless preremoulded. This strength is applied to the up-slope, initially stable, volume. We tested several s y 0 values to obtain s y 0 ¼ 4 MPa that can, in combination with the applied value ranges for the residual yield strength s y r (see Table 2) and remoulding rate C, Table 1 Default parameters  (2000) suggest an added-mass coefficient C m in the order of magnitude of 1 for short blocks. However, C m is also dependent on the ratio between landslide thickness and landslide length (Enet and Grilli 2007). C m ¼ 0:1 is used herein. For the tsunami modelling, we use a sea water density q w ¼ 1000 kg m À3 and a Manning coefficient M ¼ 0:025 . With these default paramteter values, we investigate the sensitivity to geometrical parameters such as the landslide directivity, the total and preremoulded release volume (see Table 2): 1. Landslide directivity: A constantly inclined failure plane separates the viscoplastic landslide material from the stable volcano rocks, whose dip angle of 8:2 (Giachetti et al. 2012) defines the angle of the landslide entering the water. Satellite images indicate a south-west failure direction of the 2018 event, e.g. h ¼ 135 west from the north (Babu and Kumar 2019). In order to investigate a sensitivity to the dip direction, which controls the landslide directivity, we include also scenarios with failure plane dip directions of h ¼ 125 and h ¼ 145 (see Fig. 2). 2. Total release volume: The height of the failure plane defines the total release volume, i.e. the entire viscoplastic material above the failure plane. Giachetti et al. (2012) used a total volume of 0:28 km 3 as a hypothesis for the 2018 Anak Krakatoa flank collapse. We employ this volume, as well as a smaller volume of 0:21 km 3 , to model the sensitivity to the landslide volume. Our choice of landslide volumes is consistent with the volume range of 0.22-0.30 km 3 previously reported by Grilli et al. (2019), and 0.175-0.236 km 3 reported by Heidarzadeh et al. (2020) for the 2018 Anak Krakatau flank collapse. 3. Pre-remoulded release volume: We split the total volume into a down-slope, initially weak, and an up-slope, initially stable, volume (see section above). The down-slope volume defines the preremoulded release volume that we alter from 1 to 60%, and finally 100% of the total release volume in order to indicate that a single collapse event triggered the tsunami. The volume will fail sequentially if the pre-remoulded release volume is less than 100%. Figure 3 illustrates three different values for a pre-remoulded release volume in a transect through the volcano.
Apart from the geometrical parameters, we also investigate the sensitivity to the residual yield strength (see Table 2). The yield strength is remoulded to the residual yield strength over time for materials originating from the up-slope, initially stable, volume, as defined by Eq. 1. The down-slope, initially unstable, volume is pre-remoulded, meaning that the residual yield strength is applied at the initial time. In this study, we investigate the sensitivity to the remoulded yield strength for s y r ¼ 1 kPa, 100 kPa, which are reasonable values based on preliminary tests.
In addition to the model parameters for the sensitivity studies explained above, we test the effect of friction and pressure drag coefficients, C F and C P , respectively, on the landslide velocities. We base our value ranges, C F ¼ 0:001, 0.0025, 0.005, 0.01 and C P ¼ 0:1, 0.5, 1.0, on the value ranges from studies by De Blasio et al. (2004) and Kim et al. (2019). It was found that for this particular setting, these different C F and C P values only gave differences between 4 to 5% on the maximum landslide velocities and with little influence on the tsunami generation, hence, they were not investigated further.

Applied Bathymetric Grids
Our analysis is based on the topo-bathymetric grid provided by Giachetti et al. (2012). They used data which are based on a combination of four digital elevation models, (I) ASTER topography with a 30 m-resolution, (II) submarine bathymetric maps from the Indonesian navy, (III) GEBCO bathymetry with around 900 m-resolution of the entire Sunda Strait region, and (IV) a bathymetric map of the Krakatoa Archipelago from Deplus et al. (1995). The merged topo-bathymetry has a resolution of 100 m in both x-and y-direction in Cartesian coordinates (WGS 84). The topo-bathymetry from Giachetti et al. (2012) is constructed prior to the failure. In order to produce the topo-bathymetry for our failure area without the landslide material, we insert the 8:2inclined failure plane into the volcano at various heights and orientations, as mentioned above, and cut the material above the failure plane.
In the applied bathymetry, north of the volcano Anak Krakatoa, there is a shallow water basin that has a depth of approximately 10 m. Because the simulations are highly sensitive to the accuracy of the bathymetry, and we observe significant offsets between simulated and observed dominating wave periods and arrival times, we cannot rule out inaccuracies in the applied bathymetry. We remark that this mismatch was also observed by the simulations carried out by Grilli et al. (2019) and Paris et al. (2020). We also note that such errors in open source data previously have been shown to be significant for other areas in this region, such as Palu. Therefore, we deepen the seafloor to 60 m over a total area of 174 km 2 reaching 19 km in latitudinal direction (see Fig. 1) to remove the shallow region north of the volcano. We have no firm justification for this manipulation of the bathymetry, save that it turns out to improve the agreement between observations and simulations. It may also serve as an example of the effects of uncertain bathymetry. To avoid abrupt transitions in the depth, we increase the depth between two adjacent cells in the modified area by a maximum of 3 m. We use this modified data, as well as the original data from Giachetti et al. (2012), for a sensitivity to our pre-collapse bathymetry (see Table 2). For instance, a deeper bathymetry can lead to changes in earlier arrival times, and also makes the waves less prone to dissipation and breaking in the early phase of propagation. The depth matrix for simulating the landslide dynamics covers the Krakatoa complex in a rectangle with length 15 km in longitudinal and 14:5 km in latitudinal directions, respectively. The applied grid resolution is 36 m in both directions, which is interpolated from the original topo-bathymetric data with a resolution of 100 m provided by Giachetti et al. (2012). The values of the grid cells adjacent to the boundary are copied into two ghost cells. This aids the implementation of non-reflecting outflow (LeVeque 2002). The depth matrix for the tsunami propagation computation covers the entire Sunda Strait region, including the Krakatoa complex, the western Java coast, and the southern Sumatra coast, forming a square of 180 km side length. The resolution in GeoClaw is 175 m in both directions, which is interpolated from the 100 m-resolution topo-bathymetry, and the resolution in GloBouss is 100 m. We apply non-reflecting outflow boundaries to the grid used by GeoClaw. In GloBouss we apply a sponge layer with 2:5 km at all boundaries, which relaxes the surface elevation , and to avoid artificial coastal oscillations, we apply a minimum computation depth of 1 m. The accuracy using the different model resolutions for BingClaw, GeoClaw, and GloBouss through grid refinement tests is provided in the appendix.
We analyse two regions along the western Java coast for the tsunami inundation simulations with GeoClaw. The northern region, located east of Anak Krakatoa around the Marina Jambu gauge, is defined by a rectangle with length 8:4 km in longitudinal and 11:7 km in latitudinal directions, respectively. The southern region, located south-east of Anak Krakatoa around Pantai Legon Tanjungjaya, has an extent of 7:0 km in longitudinal direction, and 16:7 km in latitudinal direction (see Fig. 4). The grid resolution in both longitudinal and latitudinal directions is 11 m for both regions, which is interpolated from the original topo-bathymetric grid with a 100 m-resolution. As the topographic data is based on open source data, the accuracy of the inundation distances can be expected to be limited. Hence, we limit our inundation analysis to compare with observed heights.

Observations
Muhari et al. (2019) studied pre-and post-event bathymetric data in the caldera formed by the 1883 eruption. Their analysis showed more than 50 m thick sediments in the 1883 caldera, which we use to validate our landslide model.
Observed tsunami surface elevation time series at the four different gauge stations, Kota Agung and Panjang at the southern coast of Sumatra, and Ciwandan and Marina Jambu at the eastern coast of Java (Muhari et al. 2019), serve as comparison to our simulated surface elevation time series. Each gauge station is equipped with three sensors. We set our initial sea level relative to Mean High Water (MHW). In Figs. 9, 10, 12, 13, and 14 we show detided surface elevation time series in Kota Agung, Panjang, and Ciwandan from sensor 2. All three sensors recorded very similar data. For Marina Jambu we show the detided data from sensor 1 and 3, because sensor 1 measured a leading crest twice as large as sensor 3, and sensor 2 did not work when the waves arrived. The sampling rate of the observed tidal data is once a minute, which is quite coarse given the 5-15 min observed first arrival wave periods (see Table 3).
Due to the coarse resolution of the shoreline in our bathymetry, we moved the gauges 300 m off the coasts. Taking into accounts the depths, that are in the range 1-6 m, say, this may shift arrival times about 100 s. This is of the same order as the deviations in arrival times that we encounter in the simulations, which illustrates the sensitivity of the arrival time to the details of the coastal bathymetry. We list measured first arrival times, surface elevations, and first wave periods of the tidal time series in Table 3. 2500 T. Zengaffinen et al. Pure Appl. Geophys.
The first arrival time is defined by the time of the first sea-level rise with a threshold value of 5 cm due to noise in the observed time series. The same extraction method is used for simulated wave metrics, but due to lower noise level, the threshold value is set to 1 cm. For the simulations the reduced threshold will make a significant difference only for the Panjang wave gauge. Here, the low 1 cm threshold leads to inclusion of the low precursor in the simulations between 3800 and 3960 s, say. This low elevation presumably has taken a different path from the source than the following crest of 0:41 m and is thus not included in the calculation of the period of the first simulated wave. While the low precursor reduces the difference in travel time between the simulation and the observation, it must be pointed out that the wave shapes are very different (see Figs. 9,10,12,13,and 14). In addition, the observed time series are coarsely resolved. Hence, the interpretations of the arrival times are uncertain. The first wave period is defined by twice the time difference between the first zero-

Assessment of Dispersion Effects
A measure of the importance of dispersion is which is a temporal variable for the evolution of dispersion effects (Glimsdal et al. 2013). Here, h is the characteristic water depth, L the propagation distance, and k ¼ c 0 T the wave length of the wave front in the source area with linear wave speed c 0 and wave period T. We evaluate c 0 and T in two gauges 10 km southwest and east of the volcano, which results in k 1 ¼ 5:16 km for the southwestern gauge and k 2 ¼ 4:31 km for the eastern gauge. We evaluate the importance of dispersion, first, for waves travelling from the southwestern gauge to Kota Agung, and second, for waves travelling from the eastern gauge to Ciwandan. Using an average characteristic water depth h 1 ¼ 320 m, propagation distance L 1 ¼ 90 km, and k 1 ¼ 5:16 km, we get s 1 ¼ 0:402 for the waves in Kota Agung, which implies a significant dispersion effect. The reason for this significant effect is mostly due to the relatively deep basin between the source and Kota Agung. As a consequence, our results for this gauge give only an indication as we do not include dispersion for the GeoClaw simulations. On the other hand, using h 2 ¼ 80 m, L 2 ¼ 45 km, and k 2 ¼ 4:31 km, the dispersion effect is small, s 2 ¼ 0:0216, for waves in Ciwandan. A similar s implies for the gauges in Marina Jambu and Panjang, thus dispersion effects are small in these shallow parts in the Sunda Strait. We run GloBouss from 600 s after the flank failure (see Sect. 2.2) with both LSW and dispersive equations. Two gauges that are both positioned two thirds of the distance from the volcano to Kota Agung and to Ciwandan (see Fig. 1), are used to investigate the importance of dispersion. Concerning the leading crest, the difference between the dispersive equations and the LSW equations is 2% and 5%, respectively, for the two gauges. For the later parts of the wave trains, the deviations are larger, in particular for the wave gauge on the way to Kota Agung (see Fig. 5). To compare the sensitivity due to the choice of model, we also compare the LSW version of GloBouss with GeoClaw at the same wave gauges. Here, the leading crest deviates by 7% between the two models, whereas later waves deviate to a larger extent (see Fig. 5), which may be due to artificial coastal reflections in the GloBouss model.
GloBouss results indicate that dispersive effects are small in shallow parts of the Sunda Strait, which is north, east, and south of Anak Krakatoa, and that a non-dispersive tsunami model is applicable. Heidarzadeh et al. (2020) reported a s in the same order of magnitude, and concluded that the effect of dispersion in the Sunda Strait is weak. On the other hand, Paris et al. (2020) found a larger influence of dispersion. However, their source model was quite different from the one used herein, and they also studied propagation in the deeper regions outside the Sunda Strait.  Figure 6 shows contour plots of the landslide masses prior to the release, 15 s, 30 s, and 90 s after the release. The landslide mass starts flowing radially over the volcano with a maximum speed of 35 m s À1 towards southwest (see Fig. 7). When the material arrives in the submarine caldera, 30 s after initiation, the landslide is guided by the shape of the caldera. After 600 s, the mean landslide speed is reduced to 2 m s À1 , and after 780 s the landslide stops moving and exhibiting a maximum deposit of ca. 50 m (not shown here). The sediment thickness is similar to the findings from Muhari et al. (2019). Figure 8 shows six colour plots of the tsunami surface elevations after 30 s, 60 s, 90 s, 120 s, 600 s, 1500 s after the landslide mass release. The waves start propagating radially inside the Krakatoa complex with a main propagation direction towards southwest. After 120 s the waves escape the Krakatoa complex through three openings between the three surrounding islands of Anak Krakatoa.

Parametric Sensitivity Study
We compare surface elevation time series for different values of the total release volume V tot , preremoulded release volume V r , landslide directivity h, residual yield strength s y r , and the bathymetry north of the volcano h b . The values of these parameters are presented in Table 2. One parameter is altered at a time, whereas the bold values from Table 2, which refers to the basic scenario, are used as default values. Also wave periods and arrival times are compared (see Table 3).

Total Release Volume
First, we show the sensitivity to the total release volume V tot to the tsunami surface elevation time series in Fig. 9 for an instantaneous collapse. No remoulding is applied here, meaning that the yield strength, s y r ¼ 1 kPa, is constant for all times and the entire collapse volume. For all four gauges, the higher volume, V tot ¼ 0:28 km 3 , causes higher surface elevations than the lower volume, V tot ¼ 0:21 km 3 . For instance, in Panjang and Kota Agung, the leading crests for the larger landslide volume are about 1.5 times higher, respectively.

Residual Yield Strength
The next parameter we investigate is the residual yield strength s y r for the instantaneous collapse with V tot ¼ 0:28 km 3 . The strength for the entire collapse volume is s y r . Figure 7 shows the time evolution of the maximum and mean landslide velocities for both s y r . Velocities for s y r ¼ 1 kPa are larger and the landslide motion lasts longer than for s y r ¼ 100 kPa. The influence of s y r on the surface elevation time series at the four gauges is shown in Fig. 10. The time series for a yield strength s y r ¼ 100 kPa show general lower surface elevations in all four gauges than for s y r ¼ 1 kPa. Overall, the smaller s y r provides better agreement with the observations, but not consistently.

Gradual Mass Release
The scenarios above had a pre-remoulded volume V r ¼ V tot , which refers to an instantaneous collapse without remoulding. We vary the pre-remoulded release volume V r \V tot (see Sect. 2.3) to simulate a gradual mass release (see Fig. 12). The down-slope volume that is pre-remoulded fails first, and the upslope, initially stable, volume is being remoulded during the landslide flow (see Fig. 11). The surface elevation time series of a gradual mass release over 240 s for V r ¼ 0:6 V tot has same wave periods and first arrivals as the instantaneous collapse, but in  Table 2 2506 T. Zengaffinen et al. Pure Appl. Geophys. general lower surface elevations, except for the Kota Agung gauge. A pre-remoulded release volume of only V r ¼ 0:01 V tot causes lower wave heights, later arrival times, and longer periods for all four gauges. The first arrival wave in Kota Agung has half the surface elevation and arrives 400 s later than the waves caused by the instantaneous collapse. In Ciwandan and Marina Jambu the first arrival waves are troughs rather than crests, and reach the gauges around 550 s later than the waves caused by the instantaneous collapse. In Panjang surface elevations are less than  Table 2 Vol. 177, (2020) Modelling 2018 Anak Krakatoa Flank Collapse 2507 0:02 m, whereas they reach up to 0:4 m for the instantaneous collapse (V r ¼ V tot ). These differences imply that the gradual mass release has a major on the tsunami generation for such a setting. In total, the instantaneous collapse yields somewhat better agreements with the observations than V r ¼ 0:6 V tot . Grilli et al. (2019), Paris et al. (2020), andHeidarzadeh et al. (2020) all based their simulations on an instantaneous collapse, which seems reasonable based on the results of our study. The very small V r in our study is definitely off the mark.

Landslide Directivity
We turn our attention back to an instantaneous mass release and focus on the landslide directivity h (measured counter-clockwise from the north). Figure 13 shows the sensitivity to h on the surface elevation time series. h ¼ 145 gives too high first waves in Ciwandan and Marina Jambu, and a too low first wave in Panjang, while there are only moderate differences between h ¼ 125 and h ¼ 135 . The first wave in Kota Agung is hardly affected by h. Later arrivals do not depend on h in a systematic manner for any gauges.

Modification of Bathymetric Depth
As a last test we investigate the water depth north of the landslide source (see Fig. 14) by applying an instantaneous volcano flank collapse. In the original bathymetry there is a shallow region with a minimum depth of h b ¼ 10 m north of the volcano (see Fig. 1). For our analysis, we increased the water depth in this region to obtain h b ¼ 60 m, which was applied in all scenarios in Figs. 9, 10, 11, 12, 13. A variation in h b has influence on the first arrivals in Panjang, whereas the first arrivals in the other three gauges are not sensitive to the change in bathymetry north of the volcano. In Panjang, the leading crest is 25% lower for h b ¼ 10 m than for h b ¼ 60 m, and the arrival time of the first steep wave front is around 100 s later for the lower h b . Subsequent crests differ in all gauges for varying h b .
Waves travel faster through the modified bathymetric area north of the volcano than through the shallower original bathymetry, which implicates an earlier first wave arrival in Panjang. Our simulated first arrival time in the Panjang gauge matches the observation better than with the original bathymetry. Moreover, waves travelling through the shallow region north of the volcano in the unmodified bathymetry are strongly steepening in the front. Analyses show that the first wave period increases with travel distance for the original bathymetry, while the first wave period is less subject to change when using the modified bathymetry. However, the simulated wave period in Panjang is still somewhat larger than the observed wave period. Also, the modified bathymetry leads to different interference patterns, which affects the waves at the Panjang gauge. Our simulated maximum surface elevation matches the observation better with the modified bathymetry (see Table 3). We find similar discrepancies in the Kota Agung gauge where our simulated first arrival wave period is larger and the leading crest smaller than the observation (see Table 3). In addition to the small discrepancies due to the significance of dispersion (see Fig. 5), we assume inaccuracies of the bathymetry in general based on the surface elevation time series mismatch between observation and simulation at the coast of Sumatra. The coarse resolution and noisy appearance of the recordings may suggest that the recordings at the coast of Sumatra are inaccurate. The lack of the complex geometries of the harbours around the gauge stations in the applied topo-  Table 2 Vol. 177, (2020) Modelling 2018 Anak Krakatoa Flank Collapse 2509 bathymetry presumably leads to mismatches (Grilli et al. 2019).

Analysis of Best Fit Scenario
An instantaneous volcanic flank collapse with a pre-remoulded and total release volume of V r ¼ V tot ¼ 0:28 km 3 (similar to the one used by Grilli et al. 2019), a landslide directivity of h ¼ 125 west from the north (similar to the one used by Grilli et al. 2019;Heidarzadeh et al. 2020), a residual yield strength of s y r ¼ 1 kPa, and a modified water depth north of the volcano of h b ¼ 60 m produces a tsunami that matches the observed surface elevation time series at the four gauge stations best. This is the basic scenario and is regarded as best in respect of a visual  Table 2 2510 T. Zengaffinen et al. Pure Appl. Geophys. comparison of wave arrival time, wave period, and surface elevation between observation and simulation. We compare our simulated wave metrics at the four gauge stations with observations (see Table 3).  Table 2 Vol. 177, (2020

Inundation Study
The best fit scenario is used for an inundation study for the two selected regions at the Java coast described above (see Fig. 4). As stated above, the topographic data used to create the grid for the inundation analysis are based on open source information. As reported elsewhere (Griffin et al. 2015), such data can often lead to hindering the onshore flow and hence limit the horizontal inundation. Hence, we do not show inundation distances. Moreover, we remark that our ability to model local inundation effects are more limited than if high resolution data were available.
The maximum simulated runup height in the northern region is 6:1 m, and 7:9 m in the southern region (see Fig. 15). In the northern region, the simulations agree well with the observations. Yet, for the three northernmost positions, the observed runup heights are underestimated by the simulations. However, the maximum simulated surface elevation in the Marina Jambu gauge is 1:4 m, which is the same as the maximum observed surface elevation. In the southern region runup heights are underestimated at all locations, except one. Here, we lack offshore tsunami data.

Concluding Remarks
We used the depth-averaged viscoplastic numerical landslide model BingClaw to simulate the 2018 Anak Krakatoa volcanic flank collapse, and coupled it to the non-linear numerical model GeoClaw that simulates the generation, propagation, and inundation of the resulting water waves. We showed that the flank collapse most likely took place as an instantaneous collapse rather than a gradual flank failure moving up-slope. For the instantaneous collapse, our simulations show a good overall match with observed wave periods, surface elevations, and arrival times at the two gauge stations on Java. On Sumatra, the simulations compare less well, and surface elevations are matched with our modified topo-bathymetric data at one gauge, and the first arrival time is matched at the other gauge. The landslide directivity is a significant control on surface elevation time series, which can be seen in three of the four gauge stations. Our best match is 125 west from the north. At the locations where we extracted runup heights from the inundation study, half the simulated data matches observed data, and the other half underestimates the observation. Shallow waters dominate the northern and eastern Sunda Strait, which makes the tsunami modelling, based on open source data, sensitive to bathymetric errors as recently observed for the tsunami in Palu Bay. According to our analysis, the first arrival time in Panjang is matched better with the modification of the bathymetry, making it deeper north of the volcano, than with the original bathymetry. However, our simulated wave periods in Panjang are still three times as large as the observed wave periods. A possible reason for the large wave periods in the model, is due to artificially high dissipation and breaking in shallow regions. This may again be due to errors in the bathymetry, but this is difficult to verify without better bathymetric data.
We acknowledge that the simulations contain many uncertainties concerning the simplified failure plane, the estimation of the total collapse volume, the depth-averaging of the slide motion, the coarse topobathymetric data, and the fact that we modelled a dacite rock debris flow with a model previously used for clay-rich sediments. The post-collapse topobathymetric data, which has recently been acquired, will serve as a base for more accurate volume estimates, a more detailed failure surface, and an extent of the landslide run-out, which can be used to validate the landslide model. That will, next to our fairly well suited model, give a more consistent description of the 2018 Anak Krakatoa volcanic collapse and tsunami.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

Appendix: Grid Refinement Tests A.1 The Landslide Model
A grid refinement test using three different spatial resolutions, Dx B ¼ 73 m, 36 m; 18 m, for the landslide model BingClaw is executed. The finest resolution Dx B ¼ 18 m is the reference resolution. The maximum landslide thickness shows a deviation of 15% between the reference resolution and Dx B ¼ 73 m, and 4% between the reference resolution and Dx B ¼ 36 m (see upper panels in Fig. 16). A resolution of Dx B ¼ 36 m is applied for further analyses. In BingClaw we employ an adaptable time step, such that the CFL numbers are always below 0.8 and with a preference value of 0.65 (LeVeque 2002).

A.2 Landslide Updating Within Tsunami Model
The source term in the tsunami model is the volume displacement from the landslide model, which is updated at a certain frequency. We tested the following update frequencies of the landslide source into the tsunami model GeoClaw: Df ¼ 12 min À1 ; 4 min À1 ; 2 min À1 ; 1 min À1 , with Df ¼ 12 min À1 being the reference frequency. The surface Vol. 177, (2020) Modelling 2018 Anak Krakatoa Flank Collapse 2513 elevations outside the source area deviate 8% between the reference frequency and Df ¼ 12 min À1 , but less than 1% for Df ¼ 4 min À1 . We apply Df ¼ 4 min À1 for further analyses (i.e. updating the landslide source terms every 15 s).

A.3 The GeoClaw Model
We tested various spatial resolutions for GeoClaw, Dx G ¼ 351 m; 175 m; 87 m, shown as time series at all four gauges in Fig. 9. The average deviation of the leading crest between the reference resolution (Dx G ¼ 87 m) and Dx G ¼ 351 m is 44% for all four gauges. The same average deviation for Dx G ¼ 175 m is 10%. The resolution Dx G ¼ 175 m represents a reasonable compromise concerning model inaccuracies and large computation times when using the reference resolution. A visual inspection of wave periods and first arrival times, in Fig. 9, suggest a relatively good match between Dx G ¼ 175 m and the reference resolution. We apply a spatial resolution of Dx G ¼ 175 m for further analyses. We refined two specific regions near the coasts (see Fig. 4) down to a grid resolution of 11 m for the inundation study. As for BingClaw, an adaptable time step is employed, this time with a CFL-number maximum of 1.0 and 0.75 as desired value.

A.4 The GloBouss Model
In GloBouss we applied grid resolutions of Dx Gl ¼ 200 m; 100 m; 50 m, with the latter being the reference resolution. For both LSW and dispersive equations, the first crest deviates with 25% between the reference resolution and Dx Gl ¼ 200 m, and with 6% between Dx Gl ¼ 100 m and the reference resolution. We apply Dx Gl ¼ 100 m for further analyses. The time step is kept constant in GloBouss. For the LSW equations we used a time step that gave a maximum CFL number of 0.60. For the disp equations the model is implicit and hence more stable. Thus, the maximum CFL number is chosen as 7.22. This gives a more favourable time step in the shallow regions of the Sunda Strait.
Publisher's Note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.