Numerical modeling of the subaerial landslide source of the 22 December 2018 Anak Krakatoa volcanic tsunami, Indonesia

The eruption of the Anak Krakatoa volcano (Indonesia) in December 2018 produced a destructive tsunami with maximum runup of 13 m killing 437 people. Since the occurrence of this rare tsunami, it has been a challenge as how to model this tsunami and to reconstruct the network of coastal observations. Here, we apply a combination of qualitative physical modeling and wavelet analyses of the tsunami as well as numerical modeling to propose a source model. Physical modeling of a volcano flank collapse showed that the initial tsunami wave mostly in- volves a pure-elevation wave. We identified initial tsunami period of 6.3 – 8.9 min through Wavelet analysis, leading to an initial tsunami dimension of 1.8 – 7.4 km. Twelve source models were numerically modelled with source dimensions of 1.5 – 4 km and initial tsunami amplitudes of 10 – 200 m. Based on the qualities of spectral and amplitude fits between observations and simulations, we constrained the tsunami source dimension and initial amplitude in the ranges of 1.5 – 2.5 km and 100 – 150 m, respectively. Our best source model involves potential energy of 7.14 � 10 13 – 1.05 � 10 14 J equivalent to an earthquake of magnitude 6.0 – 6.1. The amplitude of the final source model is consistent with the predictions obtained from published empirical equations.


Introduction
A large tsunami reaching a maximum runup height of 13 m (Muhari et al., 2019) was generated in the Sunda Strait, Indonesia due to volcanic activities of the Anak Krakatoa Volcano in the evening of 22 December 2018 ( Fig. 1). No large earthquake was reported before the tsunami whereas volcanic activities were observed at the Krakatoa Volcano from a few months earlier until the time of the destructive tsunami of 22 December, according to media reports. Therefore, it is widely believed that the tsunami was generated by volcano-related mass failures and associated subaerial landslides. Satellite images from the Sunda Strait before and after the event revealed that a large portion of the middle Krakatoa volcano island (known as Anak Krakatoa) collapsed into the sea as a result of the volcanic eruption (Fig. 1b). The tsunami left a death toll of at least 437 (https://reliefweb.int) as of 6th January 2018. The Krakatoa volcano tsunami occurred approximately three months after another destructive tsunami in east Indonesia, the 2018 Sulawesi tsunami, claiming around 2000 lives, which was partially attributed to subaerial/submarine landslides (Heidarzadeh et al., 2018;Takagi et al., 2019). The main eruption, responsible for the Anak Krakatoa tsunami, was also recorded at seismic stations in the region (Fig. 2). Based on these seismic records, the origin time of the eruption was approximately at 13:56 UTC on 22 December 2018 (Fig. 2).
Indonesia is home to approximately 130 volcanoes (Lavigne et al., 2008), a number of which are located in the sea and are capable of producing tsunamis. Possibly the predecessor of the December 2018 destructive volcanic tsunami was another deadly volcanic tsunami around the same place in 1883 generated by the Krakatoa eruption with 36000 deaths (Choi et al., 2003;Pelinovsky et al., 2005). Based on Bryant (2001), the 1883 event produced at least 20 Km 3 of pyroclastic deposits. For the 1883 event, tsunami wave runup was approximately 15 m on the average in the near field, i.e. coasts of Sunda Strait, but maximum wave runup was up to 42 m during the 1883 Krakatoa tsunami (Choi et al., 2003). Nomanbhoy and Satake (1995) modelled the 1883 event using nonlinear shallow water equations and successfully reproduced the 1883 tsunami. Latter (1981) provided a detailed review of the 1883 Krakatoa tsunami including origin time, period and coastal amplitudes of the waves and arrival times.
Tsunamis generated by volcanic activities are generally less understood due to the scarcity of such events. To reproduce the propagation and coastal amplification of volcanic tsunamis, it is essential to estimate the source of these tsunamis. Several authors in the past have studied source mechanism of volcano tsunamis. Latter (1981) listed 10 mechanisms for generation of volcanic tsunamis worldwide among which the most destructive are those generated by the impacts of pyroclastic flow deposits into the sea as well as volcanic-triggered landslides. According to Latter (1981), tsunamis from pyroclastic flows are rare while those generated by volcanic-triggered landslides are more common. Paris et al. (2014) presented a review of volcanic tsunamis including 40 incidents worldwide in the period of 1550-2007 and discussed their generation mechanisms. Eight mechanisms were proposed by Paris et al. (2014): underwater explosion, pyroclastic flow, earthquake, flank failure, caldera subsidence, air wave, lahar and collapse of the lava bench. Caldera uplift associated with volcanic activities can be considered as another generation mechanism for volcanic tsunamis as evidenced during the 2015 Torishima volcanic tsunami, Japan (Sandanbata et al.,   Tinti et al. (2006) who measured a maximum runup of 11 m and associated the tsunami to two landslides with a combined mass volume of 0.02-0.03 km 3 . Satake (2007) modelled the 1741 Oshima-Oshima volcanic tsunami in the Japan Sea whose landslide volume was estimated to be approximately 2.5 km 3 . A new modeling of the Oshima-Oshima volcanic tsunami was conducted by Ioki et al. (2019) estimating the volume of the collapsed material at 2.2 km 3 .
The December 2018 Krakatoa volcanic tsunami is among a few examples of relatively well-recorded volcanic tsunamis worldwide as the tsunami was recorded by several tide gauges. It is, therefore, an important event whose study will contribute to better understanding of coastal hazards from volcanic tsunamis. Here, we study the source of the recent Anak Krakatoa volcanic tsunami using tide gauge data and numerical simulations. We propose a static source model which fairly reproduces the observed tide gauge data.

Data and methods
Our approach to estimate the tsunami source model of the December 2018 Anak Krakatoa tsunami was based on: first, determination of the type of the source (either a dipole depression-elevation wave or a pureelevation wave) using qualitative physical experiment; second, approximation of the source dimensions using Wavelet analysis of the tide gauge records; and third, validation of the source model using forward numerical modeling of the tsunami. Our tsunami numerical modeling is based on assuming an initial tsunami wave at the end of the generation phase and then modeling its propagation and coastal amplification. In fact, we consider a static model and ignore the initial complicated tsunami generation phase. The term "static model" here implies that the complex two-phase flow consisting of water-soil interaction at the beginning of the slide motion (10-100 s) is ignored. This approach has been successfully applied to the 1998 Papua New Guinea landslide tsunami (Synolakis et al., 2002;Okal and Synolakis, 2004;Heidarzadeh and Satake, 2015a), the September 2013 landslide tsunami in the NW Indian Ocean (Heidarzadeh and Satake, 2014) and the 1945 Makran tsunami (Heidarzadeh and Satake, 2017). This forward modeling trial-and-error approach is useful for the case of the Anak Krakatoa tsunami because alternate methods, such as tsunami inversion, may not yield satisfactory results due to the complexity of the source mechanism and lack of high-resolution bathymetry data. For our trial-and-error approach, first we consider a hypothetical source and conduct tsunami propagation and compare the synthetic waveforms with actual tide gauge observations. Wee modify the assumed tsunami source according to the quality of match between synthetic and observed waveforms and continue this process until satisfactory match is achieved.
Sea level data consists of six tide gauge records (see Fig. 1a for locations) with sampling intervals of 1 min provided by the Agency for Geo-spatial Information, Indonesia (BIG) (http://tides.big.go.id). To remove tidal signals from the original records, we predicted the tidal signals using a least squares method of harmonic analysis and then removed them from the original records (Fig. 3). Wavelet, time--frequency, analysis was performed by applying the Wavelet package of   Qualitative physical experiments on the relationship between initial source length (L ini ) and tsunami wavelength (L) for two cases of volcano flank collapses at water depths of 10 cm (a) and 6 cm (b). Torrence and Compo (1998). Our Wavelet analysis was conducted using the Morlet mother function with a wavenumber of 6 and a scale width of 0.10 .
The COMCOT numerical package (Cornell Multi-grid Coupled Tsunami Model), which is based on nonlinear Shallow Water equations (Liu et al., 1998;Wang and Liu, 2006), was applied for tsunami simulations. A uniform bathymetry grid with a resolution of 8 arc-sec (approximately 250 m) was used. The bathymetry grid was interpolated from the GEBCO-2014 (The General Bathymetric Chart of the Oceans) 30 arc-sec bathymetric digital atlas (Weatherall et al., 2015). We excluded runup calculations in our numerical simulations because it requires high-resolution nearshore bathymetry/topography. Nonlinear simulations were conducted for a total time of 3 h with a time step of 0.5 s.

Type of the initial tsunami source and its dimensions
The type of the initial wave generated by volcanic tsunamis is not well understood in the existing literature. This is the reason that we had to conduct physical experiments to identify the type of the initial wave. Qualitative physical modeling was performed at the wave flume of the Earthquake Research Institute, The University of Tokyo (Japan)  to identify the type of the initial tsunami wave. This was aimed at understanding the general physics of the initial tsunami wave immediately after the generation phase rather than quantitative measurements. Froude similarity is used for physical modeling of coastal problems (Fritz et al., 2004;Sorensen, 2010); however, we do not measure wave parameters in this study because this is a qualitative assessment. The flume has dimensions of 2.6 m (length) � 0.1 m (width) � 0.3 m (height) (Sandanbata et al., 2018b). Numerous subaerial slides were made and the resulting initial waves were photographed . For submarine landslides, it was previously shown by numerous authors that the initial wave has a dipole depression-elevation shape (Synolakis, 2003;Okal and Synolakis, 2004;Okal et al., 2009;Watts et al., 2005;Enet and Grilli, 2007;McFall and Fritz, 2016). However, for the subaerial mass movements, including the flank collapse of a volcano, our experiments revealed that the initial wave is of almost pure elevation shape (Fig. 4). We conducted several experiments using the set-up shown in Figs. 4-5 by varying water depth; all of them showed initial pure-elevation wave. Two examples of our experiments are shown in Fig. 4 at two water depths of 6 and 10 cm. This was previously partially shown by Fritz et al. (2001) for the Lituya Bay subaerial slide experiments. Therefore, we consider an elevation wave as the initial wave for modeling the Anak Krakatoa Volcano tsunami. In comparison to other studies, Walder et al. (2003) also considered a wave hump (i.e. an elevation wave) to model tsunamis from subaerial landslides.
The initial estimation of the tsunami source dimension can be achieved using the equation for phase velocity of long oceanic waves (e.g. tsunamis) Rabinovich, 2010;Sorensen, 2010): in which, C is tsunami phase velocity, L is tsunami wavelength, and T is the tsunami dominating period. For long waves such as tsunamis, we know that: C ¼ ffi ffi ffi ffi ffi gh p where g is gravitational acceleration (9.81 m/s 2 ) and h is water depth (Sorensen, 2010). The tsunami wavelength (L) is α Fig. 6. Photographs and sketch showing the Anak Krakatoa volcano before and after the 22 December 2018 eruption. Photos shown in panels "a" and "b" are from https://www.jp-news.id/v/2651/dua-kali-letusan-setinggi-300-meter-tingkat-aktivitas-gunung-anak-krakatau-berstatus-waspada and www.youtube.com/watch? v¼I-3A4GR-VnU, respectively. The data of height difference for before and after the eruption is from https://sp.hazardlab.jp/know/topics/detail/2/7/27793. html. This is a Japanese website informing that the height of the mountain dropped from 338 m to 110 m. A ini and L ini are initial wave amplitude and initial length of the tsunami source, respectively.
times larger than the initial length of the source (L ini ) around the source region: Previously Heidarzadeh and Satake (2015b) reported α ¼ 2 which is confirmed by our qualitative physical experiments (Fig. 5). Therefore by combining Equations (1) and (2), we have: The sketch in Fig. 6c demonstrates the elevations of the Anak Krakatoa volcano before and after the December 2018 eruption and the initial pure-elevation wave considered in this study for modeling the propagation and coastal amplification of the tsunami. The volcano lost 228 m of its top structure to the sea due to the eruption (Fig. 6a,b,c). We consider two parameters for our static initial source model: the initial length of the source (L ini ) and the initial wave amplitude (A ini ). Twelve different scenarios were considered with varying values of L ini (1.5-4 km) and A ini (10-200 m) (Table 1) inspired by our Wavelet analyses (Section 4). A circular elevation wave is used for the shape of the initial wave (Fig. 7). For each scenario, the potential energy of the initial tsunami source is calculated using the following equation (modified Table 1 Parameters of different pure-elevation waves used for modeling initial sea level disturbance scenarios in this study.  (Fig. 6). b Initial length of the source (km) (Fig. 6). c Volume of the displaced water (km 3 ). from Satake and Kanamori, 1991): in which, E p is the potential energy of tsunami initial source in Joules, ρ is density of water which is assumed to be 1000 kg/m 3 , g is gravitational acceleration (9.81 m/s 2 ), i is a counter for the number of computational cells, n is the total number of cells in the computational domain, η i is the amplitude of the initial tsunami source scenario at the center of each cell, and S i is the area of each computational cell (S i ¼ 250 m � 250 m ¼ 62500 m 2 ). To convert the potential energy of tsunami (E p , Joules) to an equivalent earthquake magnitude (in Richter scale), the Gutenberg-Richter equation was used (Gutenberg and Richter, 1956): in which, E is energy released by the earthquake in Joules, log is the logarithm to the base 10 and M s is earthquake surface-wave magnitude.
We used E p of each source scenario as inputs for E in Equation (5). Both amplitudes and spectra of the simulated waveforms were compared with those of observations at tide gauges in order to select the best source model. For amplitude comparisons, the Normalized Root Mean Square (NRMS) misfit equation of Heidarzadeh et al. (2016a) was applied. As subaerial landslide-generated waves are usually characterized by short-period waves, it is necessary to quantify the degree of wave dispersion in order to confirm whether Shallow Water equations used by COMCOT are adequate for their simulations. We apply the following wave dispersion parameter introduced by Glimsdal et al., (2013): in which τ is dispersion parameter, λ is length of tsunami source, h is water depth around the source region, and d is the distance from source region to the coastline. Dispersion effect can be neglected for τ <0.01 but it becomes significant if τ > 0.1 (Glimsdal et al., 2013). For the case of the 22 December 2018 Anak Krakatoa volcanic tsunami, we have: λ ¼ 1.5-4 km (Table 1), d ¼ 40-80 km (Fig. 3) and h ¼ 0-100 m (Fig. 1). Therefore, we obtain an average value of 0.04 for the dispersion parameter (τ) implying that dispersion effects is not significant. Fig. 8 presents the results of Wavelet analysis which reveals that tsunami energy is distributed in the period band of 5-25 min. It can be seen that energy is channeled into varying period bands at different times. For example, three period bands of 6.3, 17 and 23 min are identified in the Wavelet plot of Marina Jambu. The high-energy channels for Kota Agung are 8.9 and 17 min. Different tsunami dominating periods are seen from one station to another. In Marina Jambu the dominating period is 6.3 min while it is 15 min in Ciwandan, 11.4 min in Panjang and 17 min in Kota Agung. Appearance of several high-energy channels and occurrence of different dominating periods can be attributed to the complicated bathymetry features in the Sunda Strait which includes shallow water (<1000 m water depth) and several small islands and coastal rias (Fig. 1). Therefore, not all of the high-energy tsunami periods can be attributed to the tsunamis source itself (e.g., Rabinovich, 1997;Heidarzadeh et al., 2016b). It has been reported by several authors that the first tsunami waves can be considered as those generated by the tsunami source (e.g., Rabinovich, 2010;Rabinovich et al., 2008;. The Wavelet analysis shows that the period of the first tsunami waves are 6.3, 7.8, 6.7 and 8.9 min in the four tide gauge stations, indicating that the average source period of the December 2018 Anak Krakatoa tsunami was 7.4 min which is the mean value of the aforesaid four periods. Applying Equation (3), by considering water depth in the range of 10-80 m and wave period of 6.3-8.9 min, we estimate the initial water surface dimension (L ini ) in the range of 1.8-7.4 km. Our source scenarios (Table 1, Fig. 7) include initial source lengths (L ini ) in the range of 1.5-4.0 km because preliminary simulations revealed that source lengths larger than 4.0 Km result is synthetic waveforms much larger than observations.

Best source model based on numerical simulations and discussion
Simulated waveforms for 12 subaerial landslide scenarios (LS) are compared with observed waveforms (Fig. 9). Results showed that scenarios with initial wave amplitude (A ini ) of 10-20 m (LS-11 & LS-12) cannot generate any waves at the location of five tide gauges. On the other hand, scenarios involving large initial length (L ini ) of 4.0 km (LS-5 & LS-6) generate much larger simulated waves than observations. In terms of quality of amplitude fit between observations and simulations, five scenarios (LS-1, LS-4, LS-7, LS-8, and LS-9) produce better performances (Figs. 9-10). These five scenarios cover initial lengths and amplitudes in the ranges of 1.5-3 km and 50-200 m, respectively. The minimum NRMS misfit belongs to LS-7 (Fig. 10).
To further narrow down our search for the best source model, we calculated the water wave spectra for the five scenarios of LS-1, LS-4, LS-7, LS-8, and LS-9 (Fig. 11). It can be seen that the spectra for the observations and simulations match fairly well in Marina Jambu and Panjang. For the two stations of Kota Agung and Ciwandan, the spectral peaks are identical (Fig. 11) although spectral gaps exist between simulations and observations. Overall, LS-9 better reproduces the observation spectra than the other scenarios.
In summary, the best model in terms of amplitude fit between observations and simulation is LS-7 whereas it is LS-9 in terms of spectral fit. Considering the uncertainties involved in this research such as lack of high-resolution bathymetry, limited tide gauge data available, and the limited number of tsunami scenarios used, it is believed that a single source model should not be chosen as the best model. Rather, the best model could have characteristics between LS-7 and LS-9 (Fig. 12). Therefore, we report the initial length and initial wave amplitude of the best model in the ranges of L ini ¼ 1.5-2.5 km, and A ini ¼ 100-150 m, respectively. Such a source model produces potential energy of 7.14 � 10 13 -1.05 � 10 14 J equivalent to an earthquake of magnitude 6.0-6.1 (in Richter scale) (Table 1). Our final source model (L ini ¼ 1.5-2.5 km; A ini ¼ 100-150 m) is consistent with published empirical equations on the maximum amplitude of the initial wave generated by subaerial landslides (Table 2) (Fritz  Noda, 1970). By using the characteristics of the December 2018 Anak Krakatoa volcanic eruption, such as average slide thickness (s) of 114 (half of the total collapsed height of 228 m, Fig. 6), average water depth (h) of 20-50 m, friction factor (f) of 0.1, slide impact angle (α) of 45 � , rock density of 2300 kg/m 3 , water density of 1000 kg/m 3 , slide volume (V s ) of 21.1 � 10 6 m 3 (volume of a cone with height of 228 m and radius of 300 m; Figs. 1 and 6), the empirical equations yield initial wave amplitude (A ini ) of 65-134 m (Table 2) which is fairly close to our A ini ¼ 100-150 m. The final initial source dimension of L ini ¼ 1.5-2.5 km is within the estimated range of source length from wavelet analysis (L ini ¼ 1.8-7.4 km). We note that the existing empirical equations for the prediction of the maximum amplitude of the initial wave generated by subaerial landslides are associated with large uncertainties as seen in the amplitude range of 65-134 m given by these equations (Table 2).

Conclusions
The 22 December 2018 tsunami generated by the eruption of the Anak Krakatoa volcano has been studied and a source model is proposed. Main findings are: � Our qualitative physical modeling of subaerial landslide failures revealed that the initial tsunami wave, immediately after the generation phase, is mostly a pure-elevation wave. Therefore, the tsunami source was assumed to be a pure-elevation wave in this study. � Based on the wavelet analyses of the tsunami observations at tide gauges, we found the period of the first tsunami signal arriving at different stations in the range of 6.3-8.9 min which guided us to an initial tsunami length of 1.8-7.4 km. � We considered 12 source scenarios, each including a pure-elevation wave, involving source dimensions of 1.5-4 km and initial amplitudes of 10-200 m. Numerical modeling of tsunami was conducted for all 12 source scenarios resulting in five scenarios as the best ones in terms of quality of waveform fits between observations and simulations. These five scenarios narrowed down the source lengths and initial amplitudes to the domains 1.5-3.0 km, and 50-200 m, respectively. Quality of spectral fit between observations and simulations further constrained the source to two models having initial lengths and amplitudes of 1.5-2.5 km and 100-150 m, respectively. � The best model (length ¼ 1.5-2.5 km; amplitude ¼ 100-150 m) involves potential energy of 7.14 � 10 13 -1.05 � 10 14 J equivalent to an earthquake of magnitude 6.0-6.1 (in Richter scale). The amplitude of our final model (100-150 m) is consistent with predictions of empirical equations which give initial amplitudes of 65-134 m for the Anak Krakatoa volcano eruption.