Time-Dependent Tsunami Source Following the 2018 Anak Krakatau Volcano Eruption Inferred from Nearby Tsunami Recordings

The eruption of the Anak Krakatau volcano, Indonesia, on 22 December 2018 induced a destructive tsunami (the Sunda Strait tsunami), which was recorded by four nearby tidal gauges. In this study we invert the tsunami records and recover the tsunami generation process. Two tsunami sources are obtained, a static one of instant initial water elevation and a time-dependent one accounting for the continuous evolution of water height. The time-dependent results are found to reproduce the tsunami recordings more satisfactorily. The complete tsunami generation process lasts approximately 9 min and features a two-stage evolution with similar intensity. Each stage lasts about 3.5 min and elevates a water volume of about 0.13 km3. The time, duration and volume of the volcano eruption in general agree with seismic records and geomorphological interpretations. We also test different sizes of the potential source region, which lead to different maximum wave height in the source area, but all the results of time-dependent tsunami sources show the robust feature of two stages of wave generation. Our results imply a time-dependent and complex process of tsunami generation during the volcano eruption.


Introduction
Large-scale volcanic events in the ocean have a great potential to trigger devastating tsunamis that could cause severe damage in nearby coastal areas. For instance, the eruption of the Krakatau volcano in 1883, one of the most explosive volcanic events in history (Maeno and Imamura, 2011), claimed over 36 000 lives and most deaths were attributed to the tsunami (Siebert et al., 2011). However, it is generally difficult to determine and simulate the tsunami generation process, due to complex phenomena associated with the volcano eruption, such as the earthquakes, landslides, and caldera collapse, etc (Latter, 1981). Previous studies have attempted to develop numerical models by categorizing the tsunamis according to the generation mechanism (Heinrich et al., 1999;Maeno et al., 2006;Imamura, 2007, 2011;Grilli et al., 2019). For example, Maeno and Imamura (2011) investigated three possible tsunami generation mechanisms following the 1883 Krakatau eruption -pyroclastic flow, caldera collapse and phreatomagmatic explosion, based on three different models -the two-layer shallow water model, the piston-like plunger model and a simple empirical model, respectively. They found that the tsunami simulations from the pyroclastic flow model matched the tsunami data well for this particular event. Generally speaking, the simulation results depend on a variety of presumed parameters, such as the flow density, duration of the eruption, and initial conditions including the geometry of the source area, the total volume and the time-dependent flux of the pyroclastic flow (Maeno and Imamura, 2011). In this study, we bypass the complexity of volcano-water interaction during the tsunami generation process, and recover the tsunami source purely based on tsunami recordings and well-developed tsunami simulation tools.
Tsunamis are mostly generated by thrust earthquakes in subduction zones, and there have been many research works on the inversion of tsunami waves to constrain the earthquake rupture process (Satake, 1987;Satake et al., 2013;An, 2015), and the associating tsunami warning strategies (An and Meng, 2016;Melgar et al., 2016;An et al., 2018aAn et al., , 2018b. If the fault geometry is relatively unclear, or the tsunamis are generated by other sources instead of earthquakes, it is more applicable to derive the initial sea surface profile (Baba et al., 2005;Saito et al., 2011;Wu and Ho, 2011). In the typical inversion method, the possible tsunami source area is divided into grids, and the initial water elevation at each grid is obtained by inverting tsunami data. Another approach to reconstruct the tsunami source is the time reversal imaging (TRI) technique, which basically adopts numerical simulations to let tsunami waves propagate from observation stations to the source area, and interfere with each other to recover the tsunami source (Hossen et al., 2015a(Hossen et al., , 2015bAn and Meng, 2017;Hossen et al., 2018). Previous studies of TRI on tsunamis assumed a static initial tsunami source, which is shown later in this study to be inferior to the time-dependent tsunami source for this event.
In this study we will present our results using the inversion approach to obtain the water elevation profile.
On 22 December 2018, the eruption of the Anak Krakatau volcano induced a destructive tsunami in the Sunda Strait, Indonesia. According to official statistics on 31 January 2019, the tsunami has caused 437 deaths, 14 059 injured and 16 are still missing (BNPB, 2019). Satellite radar images acquired ~8 h after the event reveal that the western portion of Anak Krakatau completely disappeared after the eruption, which could be a major cause of the tsunami (Williams et al., 2019). The risk of the southwest flank failure was actually warned by the previous study (Giachetti et al., 2012), yet the model underestimated the wave arrival time to Java and Sumatra coasts. Therefore, investigation of the real wave records during this event is important to improve our understanding of the mechanism of volcano-induced landslides and the associated tsunami hazards.
For the Anak Krakatau volcano tsunami, post-event field surveys have been carried out widely in Java and Sumatra islands, measuring tsunami inundation and impact (Takabatake et al., 2019;Muhari et al., 2019;Putra et al., 2020). Furthermore, the induced tsunami waves were recorded by four nearby tidal gauges, allowing for reconstructing the tsunami source process using well-developed tsunami simulation and tracing tools. Muhari et al. (2019) conducted spectral analysis of the tide gauge data to reveal the dominant periods of the tsunami,and thus estimated the source length. Heidarzadeh et al. (2020) applied a trial and error approach based on numerical simulations to constrain the tsunamis source of static initial water elevation. Grilli et al. (2019) inferred the landslide basal geometry using satellite images and aerial photography, and assumed granular material and dense viscous fluid rheologies to simulate the tsunami. They conclude that a single landslide source can explain the observed tsunami waves recorded at tidal gauges. However, seismic and infrasound data show that a separated seismic event occurred 115 s prior to the major collapse landslide process (Walter et al., 2019). Geomophological interpretation based on radar images acquired ~8 h after the event shows two discrete failure planes (Williams et al., 2019). Both studies suggest that the process of flank failure and tsunami generation is time-dependent and complex. Here, we utilize the tsunami recordings to reconstruct the temporal evolution of the tsunami source. We compare the inverted source process with seismic records and geomorphological interpretation and suggest a 9-min process divided into two main sources. The results could be helpful for further investigation of the complex volcano eruption process without knowing the source geometry and mechanism.

Data and method
In this study we adopt an inversion approach to derive the initial sea surface elevation based on tsunami recordings. The method was previously used to study the tsunami source of the 2011 Tohoku earthquake (Saito et al., 2011). In this method, the potential tsunami source region is divided into small grids, which are treated as point sources. At the center of each grid, an initial water elevation calculated from a basis function with unit amplitude is imposed. The resulting tusnami waves are numerically simulated by solving shallow water wave equations, and the waveforms at the obervation locations are obtained as Green's functions. Since nonlinearity of the tsunami waves is negligible in the propagation process, observations at tsunami stations can be expressed as linear superpositions of the Green's functions given appropriate coefficients. We adopt a least-squares method to derive these coefficients. Finally, the initial sea surface profile in the source domain is obtained from the coefficients and the basis functions.
In the 2018 Anak Krakatau event, the tsunami waves are recorded by four tidal gauges (Fig. 1). Previous studies have shown that the tsunami source caused by the volcano landslide is limited to a small region close to the failure part, and the initial sea surface dimension is about 2−8 km based on wavelet analysis and forward simulation (Muhari et al., 2019;Heidarzadeh et al., 2020). Thus,we prescribe a possible tsunami source area on the southwest of Anak Krakatau volcano, extending from 105.38°E to 105.42°E in longitude and from 6.14°S to 6.09°S in latitude. The dimension of the source area is about 4.3 km×5.6 km, divided into 6×8 grids of size 0.7 km×0.7 km. We also test two larger tsunami source areas, which are less likely to represent the actual tsunami generation region since they overlap the surrounding islands. The medium source region covers an area of about 5.6 km×6.7 km (105.38°E−105.43°E, 6.15°S− 6.09°S), and the large source area is about 13.2 km×11.2 km (105.35°E−105.47°E, 6.17°S−6.07°S). The three different source areas are shown by rectangles of different colors in Fig. 1. For each grid, we assign an initial water elevation profile using the Gaussian-shaped basis function: where x and y denote the longtitue and latitude, respectively, is the water height generated by the basis function in the ith subregion, is the center of the i-th grid, and are the width and length of the grid, and is the amplitude of the basis function we need to estimate.
The bathymetry data are extracted from the General Ba-thymetric Chart of the Oceans (GEBCO_ 2014) bathymetry with spatial resolution of 30 arc sec, which are then refined to grid size of 7.5 arc sec by spline interpolation. The linear version of tsunami simulation package Cornell Multi-grid Coupled Tsunami Model (COMCOT) (Liu et al., 1998;Wang and Liu, 2007;An et al., 2014) is adopted to simulate the tsunami propagation and calculate the Green's functions at the four tidal gauges. The total simulation time is 6 000 s, and the time step is 0.5 s to satisfy the Courant−Friedrichs− Lewy (C.F.L.) condition.

Static tsunami source
We first assume a static profile of water elevation as the tsunami source. It is known that the volcano eruption time is around 14:00 (UTC time), but the exact initial time of the tsunami source is unclear. Therefore, we conduct a search from 13:50 to 14:05 with an interval of 1 min to find the optimum initial source time. For each initial time, we adopt the non-negative least squares method in the inversion (Lawson and Hanson, 1995), on the basis that the volcano-induced landslide is assumed to cause positive sea surface displacement due to mass entering the ocean. For the three different source sizes, results consistently show that the initial time of 13:54 best matches the tsunami recordings, which is similar to the origin time of a low-frequency earthquake on the regional seismic network (  source. From Fig. 2, it is found that the initial water elevation is located to the west of Mount Anak Krakatau, which is consistent with satellite images interpretations (Williams et al., 2019). The total volumes of elevated water of the small, medium and large source sizes are estimated to be 0.03, 0.05 and 0.17 km 3 . Nevertheless, it is also found that the results of the optimum static tsunami sources are not satisfactory. The large source has the best tsunami predictions, but the initial water elevation extends to the surrounding islands, which is unlikely in reality. While the small tsunami source leads to poor tsunami waveform fit. Seismic waveforms and infrasound records both indicate that a separated high-frequency event ~115 s prior to the main event, which was interpreted as the seismic precursor or even trigger of the main sector collapse by Walter et al. (2019). Also, the spectrogram analysis indicates that the collapse process includes a 1-2-minute-long low-frequency signal followed bỹ 5 minutes of strong emissions (Walter et al., 2019). Thus, it is inferred that the time-dependent evolution of the tsunami source could be non-negligible, which can be taken into account to improve the tsunami data fit.

Time-dependent tsunami source
To account for the time-dependent evolution of the tsunami source, the possible duration time of the tsunami source between 13:50 and 14:05 is discretized to time segments with an interval of 0.5 min. For each time segment, we assume an independent tsunami source. Wave propagation from each tsunami source can be linearly added up to obtain the predicted waves at the tidal gauges. The Green's functions for all the time segments are the same except different time delays. Again, we use the non-negative leastsquares method to constrain all the tsunami sources varying in time. For a given time, the evolution of all the tsunami sources before this moment gives the water surface profile at this moment.
The evolution of the sea surface profile between 13:50 and 14:05 is plotted in Fig. 3. Note that here we give the results of the small source area, which is more likely to represent the actual tsunami source. We select only a few key moments in the evolution process displayed in Fig. 3, to highlight the two-stage feature of the tsunami source. A more detailed evolution process is provided by Fig. S7 in the Online Resource with uniform time interval of 0.5 min. Fig. 3 shows two clearly separated stages of water elevation, which appear approximately at 13:53 and 13:59, respectively. The location and height of the water elevation are similar for the two stages. The sea surface to the southwest of Mount Anak Krakatau is significantly elevated from about 13:52 to 13:55, with water height of about 3 m. The sea surface then rises again approximately from 13:58 to 14:00, with the height about 6 m. It should be noted again that the sea surface profile in each plot is the summation of the propagation of all the previous tsunami sources, instead of the newly-generated tsunami source at this moment. However, since the two stages are clearly separated, it can be inferred that the second stage is not the propagation effect of the first stage.
We also calculate the volume of newly-generated water elevation, as shown in Fig. 4. From Fig. 4, it is seen that the volume of the water elevation generated is negligible before 13:52 and after 14:00. Thus, the tsunami generation process lasts about 9 min, approximately from 13:52 to 14:00. Additionally, it is separated to two stages at about 13:56. The first stage lasts about 4 min from 13:52 to 13:55 and the second one from 13:58 to 14:00. The volume of wa- ter elevation in the two stages is estimated to be 0.13 km 3 . The volume of elevated water in time using medium and large tsunami source areas are also provided in Fig. 4 for comparison. It shows that, although the spatial distribution of the elevated water varies due to different prescribed source areas, the volume of the elevated water and its temporal evolution present high consistency. The feature of two-stage water elevation during the tsunami generation is observed regardless of the prescribed source area.
In Fig. 5, the predicted tsunami waves at the four tidal gauges are compared with the tsunami recordings, as well as the predictions from the optimum static tsunami source. At stations Panjang and Kota, the arrival time is better predicted by the time-dependent source. At stations Kota, Serang and Ciwandan, the time-dependent source predicts higher first-wave height. Particularly at Station Serang, the time-dependent source also recovers the first trough to some extent; the second crest and some of the trailing waves, which are not included in the inversion, are also better predicted. The waveform at Panjang Station is poorly predicted by both static and time-dependent tsunami sources. Previous studies also yield poor simulated results (e.g Grilli  et al., Heidarzadeh et al., 2020). We infer that the tide gauge panjang is surrounded by structures like a marina or a harbour, and these strucures are not resolved by the coarse bathymatry data adopted in this study, leading to inaccurate simulations. In summary, the time-dependent tsunami source produces improved waveform fitting than that of the optimum static tsunami source.

Discussion
Another approach of utilizing tsunami data to reconstruct the tsunami source in the absence of source mechanism is the time reversal imaging (TRI) technique. For this event, we have also attempted to apply the TRI method. The time reversal images from the four stations are first obtained separately (Panels (a) to (d) in Fig. 6), and then combined linearly to produce the final profile of initial water elevation (Panels (e) and (f) in Fig. 6). However, we find that the tsunami waves from the four stations do not interfere in the source area, and hence a satisfactory initial profile is not obtained. This is mainly because of two reasons. First, the simulation of the tsunami waves is highly affected by the local bathymetry. Among all the four stations, only the leading waves from Station Serang are clear in the source area. The waves from the other three stations are contaminated by the coastlines. Second, as shown in the previous section, the generation of the tsunami source lasts about 10 min, which is comparable to the wave period, so it is not suitable to neglect the time-dependent process of the tsunami source and recover a static image of initial water elevation. Thus, the TRI results of the remaining stations do not interfere in the source area to construct a focused static tsunami source.
In this study, we have excluded wave dispersion and solved the linear shallow water wave equations in the numerical calculations. Since the size of the tsunami source is relatively small compared with tsunamis generated by large subduction earthquakes, it is questionable if the wave dispersion is indeed negligible. Here we use the optimum static tsunami source to simulate the tsunami waves, and compare the numerical results with and without dispersion. The simulated waves at the four stations are plotted in Fig. 7. The numerical simulation with wave dispersion is carried out using software package FUNWAVE (Shi et al., 2012;Kirby et al., 2013). Although there are some high-frequency oscillations that only exist in the non-dispersive simulation, the overall simulated waves are very similar. A possible reason is that the water depth in the source area is relatively shallow. The effect of wave dispersion depends on the ratio of wavelength and water depth, i.e., kh, where k is wavelength and h is water depth. Thus water waves in shallow water generally have insignificant dispersion. Another possible reason is that the computational domain in this study is small and the four tidal gauges are near the source. The effect of wave dispersion cumulates with propagation distance, so the cumulated dispersion is negligible in short distance.
The SAR image acquired ~8 h after the tsunami reveals two discrete failure planes: one indicates the failure of the western flank, and the other indicates a surface break close to the preexisting crater (Failure planes A and B, respectively, Williams et al., 2019). Seismic and infrasound re- cords show a separated seismic event followed by a Mw 5.3 event of longer duration (Walter et al., 2019). Based on the results of time-dependent tsunami inversion, the first seismic event could correspond to the slide on the listric fault plane near the preexisting crater (Failure plane B), leading to the first stage of tsunami generation. Although the seismic derived origin time (Walter et al., 2019) is about 2 min later than the first generation of the tsunami (~13:52, Fig. 4), it stays between the two water elevation stages, and should be within the uncertainties of both inversions. The focal mechanism of the second event demonstrates that it occurs on a south-west dipping fault plane with opening mechanism, which probably represents the loss of the western flank (Failure plane A), leading to the second stage of tsunami generation. The origin time of the second event obtained by Walter et al. (2019) (around 13:56) is consistent with the start of the second tsunami generation in Fig. 4. Seismic records suggest that the duration of the second seismic event is longer than that of the first event, and it also radiates more low-frequency seismic waves (Walter et al., 2019), while our results indicate similar duration and intensity of two stages. Thus, it is still not fully understood how the first seismic event generates similar tsunami waves as the second event. Nevertheless, our tsunami wave inversion together with seismic and geomorphology analysis is in favor of a time-dependent and complex tsunami source rather than a single pulse-like source (Grilli et al., 2019).

Conclusion
In this study we use the tsunami recordings at four tidal gauges to reconstruct the tsunami source following the 2018 Anak Krakatau eruption. If assuming a static initial profile of water elevation, the optimum tsunami source occurs at 13:54 (UTC time), and the volume of water elevation is about 0.03 km 3 . However, we find that the static tsunami source does not satisfactorily predict the tsunami waves. By accounting for the time-dependent evolution of the tsunami source, results show that the tsunami generation process lasts about 8 min, from 13:52 to 14:00. It is divided to two stages, the first lasting approximately from 13:52 to 13:55 and the second from 13:58 to 14:00. The volume of water elevation in the two stages is estimated to be 0.13 km 3 for each stage. We have also tested different potential source areas in the inversion, and the results of the water volume and its temporal evolution are similar. We note that these findings are obtained solely based on the tsunami recordings. The two-stage source process in general agrees with seismic and geomophological analysis. Our results can be also integrated with more geological, seismic and geomorphological evidences if one is to extend the results to interpret the 2018 Anak Krakatau eruption process.