Modelling of the tsunami from the December 22, 2018 lateral collapse of Anak Krakatau volcano in the Sunda Straits, Indonesia

On Dec. 22, 2018, at approximately 20:55–57 local time, Anak Krakatau volcano, located in the Sunda Straits of Indonesia, experienced a major lateral collapse during a period of eruptive activity that began in June. The collapse discharged volcaniclastic material into the 250 m deep caldera southwest of the volcano, which generated a tsunami with runups of up to 13 m on the adjacent coasts of Sumatra and Java. The tsunami caused at least 437 fatalities, the greatest number from a volcanically-induced tsunami since the catastrophic explosive eruption of Krakatau in 1883 and the sector collapse of Ritter Island in 1888. For the first time in over 100 years, the 2018 Anak Krakatau event provides an opportunity to study a major volcanically-generated tsunami that caused widespread loss of life and significant damage. Here, we present numerical simulations of the tsunami, with state-of the-art numerical models, based on a combined landslide-source and bathymetric dataset. We constrain the geometry and magnitude of the landslide source through analyses of pre- and post-event satellite images and aerial photography, which demonstrate that the primary landslide scar bisected the Anak Krakatau volcano, cutting behind the central vent and removing 50% of its subaerial extent. Estimated submarine collapse geometries result in a primary landslide volume range of 0.22–0.30 km3, which is used to initialize a tsunami generation and propagation model with two different landslide rheologies (granular and fluid). Observations of a single tsunami, with no subsequent waves, are consistent with our interpretation of landslide failure in a rapid, single phase of movement rather than a more piecemeal process, generating a tsunami which reached nearby coastlines within ~30 minutes. Both modelled rheologies successfully reproduce observed tsunami characteristics from post-event field survey results, tide gauge records, and eyewitness reports, suggesting our estimated landslide volume range is appropriate. This event highlights the significant hazard posed by relatively small-scale lateral volcanic collapses, which can occur en-masse, without any precursory signals, and are an efficient and unpredictable tsunami source. Our successful simulations demonstrate that current numerical models can accurately forecast tsunami hazards from these events. In cases such as Anak Krakatau’s, the absence of precursory warning signals together with the short travel time following tsunami initiation present a major challenge for mitigating tsunami coastal impact.

to have been accompanied by an increase in summit elevation. The map in Abdurachman et al. (2018) suggests a summit height of 310 m. A height in excess of 300 m is also suggested by media reports following the collapse (e.g., 338 m referred to in 2 ), although the basis of these measurements is unclear. Photographs from the summer of 2018 3 show that the active vent was south-west of the cratered summit of the 2014 cone, and had formed a rounded pyroclastic cone on this side that projected well above the crater rim of the older cone. Based on this, we have estimated that the summit of Anak Krakatau had a maximum elevation of 330 m prior to its collapse, but note that this is not precisely constrained, and may have been closer to 300 m.
For the tsunami modelling, we used a 100 m resolution grid based on Giachetti et al. (2012), but increased the elevation of the Anak Krakatau pyroclastic cone to 330 m (by multiplying elevations by a factor of 1.1), as indicated in Fig. 3. The relatively low resolution of the subaerial island in this grid gives a smoothed coastline and edifice shape relative to the ASTER/SRTM digital elevation models, but there is little difference in the volume of the edifice or of the subaerial landslide mass. Given other uncertainties in the landslide basal plane and total volume (discussed below), we suggest that the simplification in the shape of the subaerial volcano is not likely to be a significant source of error in the tsunami source model, particularly compared to uncertainty in the pre-collapse geometry of the submarine slopes of the volcano.

S1.2. The subaerial landslide headwall
Determining the shape of the subaerial lateral-collapse headwall is key to obtaining an accurate estimate of the primary landslide volume, and particularly how much of the 330-m high pyroclastic cone was involved in the collapse. In addition to the coastline shape observed in Sentinel-1A and ALOS-2 SAR images from the 23 rd and 24 th December (Fig.  2, S.1.), aerial photographs taken on 23 rd December provide important constraints and the northern and southern margins of the collapse scar. Images posted on several websites of Surtseyan explosions on 23 rd December also show the remains of the former coastline 4,5,6,7 . This indented coastline was the result of multiple overlapping lava-delta lobes, and the margins of the scar are easily located relative to this complex coastal shape. The shape of the headwall interpreted from these images is shown in Fig. 2 and is the basis of our estimated landslide geometry. The shape has been simplified in the smoothed 100-m grid used for our simulations, while maintaining the same overall dimensions.
The Sentinel-1A image from 23 rd December shows a relatively linear, wide-angled collapse scar, but it is ambiguous whether this entirely encompasses the pre-collapse pyroclastic cone of Anak Krakatau. The bright areas in the central part of the collapse scar suggest that parts of the NE flank of the cone may remain in situ, while the southernmost part of the 1930s tuff cone crater also appears to be present. Thus, although the scar clearly cut beyond the active vent site and removed the majority of the cone, a residual part may have remained in place. This remaining material is likely to have been highly unstable, and subsequent images suggest erosion of the collapse scar to the east and north-east, cutting back towards and possibly beyond the line of the 1930s tuff cone crater (see interpretation in Fig. 2). Note however that interpretation of this early radar image is complicated by the significant radar reflectivity of the eruption plume, which partly obscures the eastern part of the island, and also the potential for radar signals to be backscattered by pyroclastic material that may have accumulated around the vent site (the potential presence of floating debris is suggested by aerial photographs taken on 23 rd December, noted above 4,5,6,7 ). Nevertheless this first radar image is the most useful one, because later satellite images also highlight the very rapid post-collapse infilling of the primary collapse scar, and growth of the island coastline to the east (Fig. 2). As a result, satellite images from after 24 th December are a poor guide to the initial collapse geometry.

S.1.3. Constraints on the landslide basal surface
The observation of Surtseyan explosions the day after the lateral collapse indicates that the vent site was submarine, and that the headwall slope was thus relatively steep (cutting from ~150 m (the estimated height of the post-collapse scar) to beneath sea level, over a lateral distance of ~300 m). As discussed in the main text, the relationship between Surtseyan explosion styles and water depth is not well characterized, but historical observations of similar behaviour occurred at vent depths ≤50 m (Moore, 1985;Belousov and Belousova, 2000;Cole et al., 2001). The location of the post-collapse crater lake, which formed as the island re-grew through to early January 2019, coincides with that of the pre-collapse vent: we thus infer that the same conduit remained active throughout the event. From this, we interpret that the collapse cut the active conduit at a depth of < 50 m but beneath sea level, and have adopted a nominal depth of 25 m in our estimated basal surface geometry. From this, it is clear that there must have been a submarine component of the landslide, but there is limited evidence on which to estimate the extent of this.
The shape of the submarine edifice of Anak Krakatau is evident from bathymetry in Deplus et al. (1995), showing that the SW flank of the cone overlaps with the margin of the 1883 caldera, forming a topographic bulge on the caldera rim. We have assumed that the initial failure of the landslide mass did not extend beyond or below the base of the Anak Krakatau edifice, and have thus selected a shape that continues from the subaerial scar, but remains within the contours of the SW submarine flank to define a smooth symmetrical form. The shape derived from this is shown in Fig. 3. If the foot of the landslide is to remain within these contours, it is clear that the steep headwall gradient must transition into a much more gentle basal surface, defining a concave basal landslide geometry that is typical of volcanic lateral collapses in general (e.g., Glicken, 1996;Watt et al., 2014). The precise shape of this is poorly defined. A more conservative, shallower basal plane would result in a reduced landslide volume (a reduction of 20% is indicated in Fig. 3).
Although bathymetric information from Anak Krakatau is fairly limited, the growth of the island has been repeatedly surveyed since 1930, and Sudradjat (1982) provides a cross section that indicates this growth pattern. Surfaces from this cross section are shown in Fig. 3 and demonstrate that our estimated landslide surface falls within the Anak Krakatau edifice, rather than cutting into older caldera-wall rocks. The surfaces also suggest that the landslide headwall did not cut back as far as the northeast rim of the pre-1960 tuff cone crater, which increased in elevation until the 1960s, when the Anak Krakatau vent became subaerial and phreatomagmatic activity ceased.

S.1.4. Uncertainties and limitations
Although the information summarised above constrains the shape of the initial landslide headwall and also highlights post-collapse modifications of the subaerial island, it is important to recognise limitation in these constraints. The shape of the central part of the headwall, and particularly its gradient, cannot be precisely constrained. The Surtseyan post-collapse activity shows that there must have been a submarine component to the landslide, but this same activity rapidly obscured and infilled parts of the scar, and may also have been accompanied by erosion of the steeper parts of the headwall. This postcollapse modification means that the high-resolution optical images taken from December 29 onwards do not help constrain the primary landslide scar geometry. To better estimate the submarine surface of the landslide is likely to require geophysical surveying.
The available images also do not allow us to determine if the collapse occurred in a single stage or in a piecemeal manner. However, both wave observations and seismic records suggest a rapid single-stage collapse, or very closely spaced stages, such as occurred during the collapse of Mount St Helens in 1980 (Glicken, 1996), and as inferred from the single observed tsunami wave train generated by the Ritter Island collapse, in 1888 (Ward and Day, 2003;Day et al., 2015). Our model results also imply that only a collapse of this type would be capable of generating the observed magnitude of tsunami. Additionally, our observations can't discern if there was an explosive or eruptionassociated component to tsunami generation, in addition to the mass movement. However, we emphasise that our model results currently suggest that such a process is not required as source component of the tsunami, given uncertainties in the collapse volume. Note that, besides using NHWAVE and FUNWAVE-TVD, some preliminary simulations and sensitivity analyses to the AK collapse parameters were performed using TSUNAMI-SQUARES, a simplified but more computationally efficient model simulating submarine/subaerial slides and edifice collapse, based on momentum exchange between discrete elements (e.g., Xiao et al., 2015).

a. b.
c. d.  Fig. 3). Reference level is MSL + 1.5 m. Maps in a. to d. were produced by the authors based on the bathymetry/topography dataset, using MATLAB version 2017b; topography from GOOGLE EARTH georeferenced satellite images was embedded in these maps using an API key: GOOGLE EARTH, 2019 DigitalGlobe.
Finally, 3D animations of NHWAVE simulation results, for viscous and granular slide scenarios are shown in files: Krakatau_slide3D_gran.mp4 for the granular slide; (ii) Krakatau_slide3D_visc.mp4 for the viscous slide; (iii) Krakatau_wave_slide3D_gran.mp4 for the granular slide with wave generation; and (iv) Krakatau_wave_slide3D_visc.mp4 for the viscous slide with wave generation.
These animations clearly show where the slide material runs out onto the seafloor, as well as near-field wave generation, propagation and interaction with the surrounding islands. Each of the tide gauges 1-4 used in this study (Ina-COAP, 2019;Fig. 1a; Table 1) is located near or within harbours/marinas or close to/behind a coastal structure (Fig S.4). Tide gauge 1, for instance, where the largest waves were recorded (Fig. S.4a), is located inside and at the far end of Marina Jambu (Fig. S.4a), which is about 200 by 100 m and not modelled in Grid 1 (especially at a 100 m resolution), within which the incoming tsunami would have triggered seiching oscillations. Assuming a quarter wavelength resonance within the shallow marina, an average depth of 3 m, the resonance period would be 2.5 min. Comparing the modelled and observed time series at gauge 1, additional sustained higher frequency oscillations (and hence less decay) can be seen in observations, particularly later in the measured data, than in the simulated data.  Table 1, Figs. 1a, S.4) at 1 min intervals, using one of the functioning sensors (blue symbols).
The red line shows the data filtered using a Butterworth filter (0.025, 0.00694 normalized passband/stopband frequencies, and 0.04 and 40 dB of passband/stopband attenuation). Time t = 0 is the start of the event estimated at 20:57' local time (UTC + 7h) on December 22, 2018.
The raw surface elevation data, recorded at each tide gauge at 1 min intervals using 2 functioning sensors (Ina-COAP, 2019), was detided using a Butterworth filter and the Matlab filtfilt function. Based on the sampling period of 1 minute for the gauge data, the Butterworth filter parameters were set assuming 20 min and 12 h tsunami and tide periods, respectively, passband and stopband normalized frequency cutoffs of 0.025 and 0.00694, and 0.04 and 40 dB of passband and stopband attenuation, respectively. Note, the normalize frequencies vary in [0, 1] rad/sample. The raw and filtered data are shown in Fig. S.5. The measured tsunami elevations shown in Fig. 5 were obtained by subtracting the filtered data from the raw data. Although simulations were limited to about 130 min (Fig. 5), Fig. S.5 shows that significant waves were recorded for a much longer time at most tide gauges, due to multiple wave reflections and re-reflections in the Sunda Straits.   Table 1) for each of the tested collapse rheology (viscous or granular) and for 3 different collapse volumes, the likeliest 0.27 km 3 , a larger 0.30 km 3 (+10%), and a smaller 0.22 km 3 (-20%) volumes (see Fig. 3). While larger differences between various cases occur in the near-field, at wave gauges located near AK, only small differences can be seen between the 6 cases in the far-field, at tide gauges, particularly early in the time series.   In the simulations the north and northwest coasts of Panjang experienced 10-18 m runups, while the southwest and southern coasts experienced 10-23 m runups (Fig. S.7.b). These predictions appear to be supported by Figs. S.8a,b, although this will need to be confirmed by actual field surveys. In these figures, the northwest tip of the island shows signs of erosions and of some overtopping flow, which was strong enough to uproot most of the vegetation on the cliff-like coastline but not atop of it. There are also widespread signs of ash deposition. This is consistent with the more moderate tsunami generation on this side of AK (see animations in Suppl. #2).
In the simulations the northern coast of Rakata experienced 30 to nearly 50 m runups, while the western coast experienced 15-30 m runups ( All of these observations will ned to be confirmed and quantified in actual field surveys, which the authors are planning to conduct in the near future.