Seismic signature detection during the 2018 Anak Krakatau flank collapse and tsunami using seismic amplitudes from regional-scale monitoring

On 22 December 2018, the major flank collapse of Anak Krakatau volcano generated a tsunami that struck the surrounding coasts of Java and Sumatra islands in Indonesia without warning. It was later suggested that the corresponding seismic event lacked the body-wave arrivals typical of tectonic earthquakes, causing difficulties for the automated detection system to recognize the event. We explore the possibility of detecting the seismic signature of such events without relying on the arrival times of body waves, by measuring seismic amplitudes in a regional seismic network at the expected arrival times from a fixed, potential source and comparing them to the theoretical attenuation of surface waves. We propose a fast detection method and evaluate the method using seismograms recorded during the flank collapse and tsunami episode as well as several known tectonic earthquakes. Detailed examinations of the detection results confirm the seismic signatures of the flank collapse and teleseismic events as suggested by previous studies. We also find more seismic signatures suggesting the occurrence of two possible smaller collapse events and variations in the eruptive activity related to the major flank collapse, as well as body and surface wave signals from two teleseismic earthquakes that were present during this episode. Finally, we construct a timeline of events during this devastating episode, combining our results with previous studies as well as insights from weather radar observations. With the ability to detect and discriminate various types of seismic events from each other, the detection method can be useful in assisting the existing monitoring and early warning systems in detecting major volcano-related tsunamigenic events.


Introduction
The devastating 22 December 2018 flank collapse and tsunami of Anak Krakatau volcano in Indonesia was one of the most significant volcano-related disasters in modern time and the latest example of tsunami generation by partial collapse of a volcanic edifice.There were many lessons learned, particularly on the application of local tsunami early warning systems, which are usually optimized to detect tsunamis caused by tectonic earthquakes.Public tsunami warnings in Indonesia are released by the Agency for Meteorology, Climatology, and Geophysics (BMKG), which typically raises an alert when a potential tsunami-generating earthquake is detected (Murjaya et al. 2012).However, despite the knowledge that Anak Krakatau was in a state of unrest, there was no alert, even after the tsunami arrived along the coastlines of western Java and southern Sumatra.The generating mechanism of the tsunami was not immediately known, since there was no potential earthquake detected around the impacted areas.Later analyses revealed the existence of a seismic event associated with the flank collapse with a magnitude ranging from M w 4.8 (BMKG earthquake catalog, https:// repog empa.bmkg.go.id/) to M w 5.9 (Ye et al. 2020) that originated in the vicinity of Anak Krakatau at around 20:55:48 WIB (Western Indonesian Time, UTC + 7), or about 30 min before the earliest tsunami arrival in western Java.Walter et al. (2019) and Ye et al. (2020) suggested that the higher-frequency body waves (P and S waves) of this event are unusually small in contrast to a tectonic earthquake with a similar magnitude.This may explain the absence of initial seismic detection of the flank collapse event, resulting in no alert being issued to the public.
Various approaches have been developed in seismic monitoring, detection, and location determination, motivated by the difficulties of observing the arrival onsets of body waves in some types of seismic events using conventional arrival time-based detection routines.These approaches provide alternatives for detecting and locating these "difficult" events by also addressing several signal-related issues such as low signal-to-noise ratio (SNR), emergent onset, and irregular variations in amplitude and duration, such as in the case of tremor-like signals.One approach is to use coherent seismic signals over many seismic stations.Array analysis (e.g., Rost and Thomas 2002) uses the plane wave assumption of distant earthquakes to measure the apparent velocity and incoming direction of seismic waves at a dense local seismic network.For an epicenter located within the network, the sum of seismogram amplitudes at the expected arrival times from trial sources (Kao andShan 2004, 2007) or a further step using the cross-correlation of seismogram envelopes (Obara 2002) have been used to locate tectonic tremors.The rapid developments of these alternative approaches have been in the field of volcano monitoring using local seismic networks, by employing seismogram cross-correlation to locate seismic sources (Droznin et al. 2015;Permana et al. 2020) or further computing the seismic covariance matrix to make use of the extracted eigenvalues (Seydoux et al. 2016) and eigenvectors (Soubestre et al. 2019;Zhu et al. 2021) for detection and location.Another approach is to locate seismic sources using the spatial distribution of seismic amplitudes, comparing them to the theoretical ones computed from a trial source (Battaglia and Aki 2003;Ogiso et al. 2015) and at the expected arrival times (Kumagai et al. 2009(Kumagai et al. , 2010) ) following a known seismic attenuation model.Some studies have taken a further step by combining the cross-correlation and amplitude spatial distribution (Maeda and Obara 2009;Permana and Aoyama 2023).
In this study, we follow these alternative approaches to investigate whether the Anak Krakatau flank collapse and tsunami episode should be detectable in near real-time or not, in the case of no clear detection of body wave arrivals at the surrounding regional seismic network.We take inspiration from the studies applying these approaches using regional networks to monitor low-frequency tremors in Japan (Maeda and Obara 2009;Obara 2002) and a mass flow event in India (Cook et al. 2021).Another example from the field of volcano monitoring is volcanic tremor studies in Kamchatka (Droznin et al. 2015;Journeau et al. 2022), while most other volcano monitoring studies are usually conducted using local networks.Because the scale of our regional network is much larger compared with the scale of the seismic source changes at Anak Krakatau (e.g., source migration, etc.), we do not focus on detailed source localization.Instead, we fix the source at Anak Krakatau and only focus on whether a surface seismic event originated at the volcano or not.We develop, evaluate, and propose a detection method based on the alternative approaches using the spatial distribution of seismic amplitudes.Our detection method is useful to monitor a known potential surface seismic source, such as a volcano, using a regional monitoring network.This is particularly useful when monitoring using local networks is difficult to perform.The method may also be useful for improving real-time monitoring and early warning developments for detecting volcanic tsunamigenic events, such as in the case of Anak Krakatau.In the following descriptions, all time information is presented in WIB (UTC + 7).

Observation data
Anak Krakatau, located in Sunda Strait between Sumatra and Java Islands (Fig. 1a), is the youngest volcanic cone that was formed after the major eruption and calderaforming collapse of Krakatau volcano in 1883 (Self and Rampino 1981).The young volcanic island is very active and dynamically changing, a condition that may pose a risk to local monitoring instruments in addition to its isolated location.In the case of major events such as the 2018 flank collapse, seismic signals may travel further and be recorded by the surrounding regional seismic stations used for earthquake monitoring.Those stations are a part of Indonesia's nationwide broadband seismic network maintained by BMKG, from which we select a subset of 13 stations that were available during the flank collapse event with epicentral distances of < 400 km from Anak Krakatau.These stations are distributed in the southern part of Sumatra, Enggano Island, and the western part of Java, with the closest station, CGJI, located approximately 64 km from the volcano (Fig. 1a).We collect the vertical ground velocity data from 20:30 to 22:30 on 22 December 2018, capturing the time of the flank collapse and the following tsunami arrivals (Fig. 1b).For the purpose of validating the detection method, we also obtain vertical seismograms of four tectonic earthquakes located by BMKG in 2018 within the Sunda Strait, south of western Java, and west of Sumatra (Fig. 1a and Table 1).We correct all seismograms for instrument response, taper both ends of the seismograms using Tukey window, equalize the sampling frequency to 20 Hz, and apply a 0.02-1 Hz bandpass filter.The seismograms during the flank collapse and tsunami episode show several visually distinct events (Fig. 1b).The first event observed around 20:55 has been associated with the flank collapse that generated the tsunami (e.g., Walter et al. 2019;Ye et al. 2020).The second event recorded around 21:35 has been suggested as teleseismic signals from an M w 6.0 earthquake in Vanuatu, followed by low-frequency signals interpreted as the tsunami arriving at the coast (Perttu et al. 2020).
We look for insights on surface phenomena during the flank collapse and tsunami episode that are difficult to observe due to the isolated location and the nighttime occurrence.Satellite-based observations have been used in previous studies to determine the cause and evolution of the flank collapse (e.g., Perttu et al. 2020;Walter et al. 2019;Williams et al. 2019).We take into consideration the direct observations based on eyewitness accounts from the surrounding waters and inland areas from Perttu et al. (2020), although the view of Anak Krakatau was reported to be obscured around this episode.Meanwhile, volcanic plumes during eruptions can also be observed and characterized through Doppler radar observations (e.g., Marzano et al. 2006Marzano et al. , 2013)), where their internal structure can be mapped in high detail (Maki et al. 2021).Radar observation of the eruption plume remains the least discussed, as there is no dedicated study at Anak Krakatau and there have been only a few applications of this monitoring approach in Indonesia (e.g., Syarifuddin et al. 2019).
Fortunately, Anak Krakatau is within the scanning range of a C-band single-polarization Doppler weather radar operated by BMKG in western Java (Fig. 1a).Therefore, we can explore this potential approach in terms of detecting volcanic eruptions through regional-scale monitoring.The radar data consist of radar reflectivity factor values (in dBZ unit) in conical scans with a 200 km radius at nine elevation angles from 0.5° to 19.5°,  preprocessed and gridded with a 500 m resolution (Permana et al. 2019).We only use elevation angles of > 0.5° to avoid the effect of electromagnetic interference at the lowest angle.The minimum temporal resolution is 8 min between consecutive observations, which is higher than satellite-based observations in the previous studies.To complement the radar data, we collect horizontal wind velocity data from the ERA5 reanalysis dataset provided by the European Centre for Medium-Range Weather Forecasts.The dataset is provided globally on a 0.25° grid at various pressure levels with a 1-h temporal resolution (Hersbach et al. 2020).We select four grid points closest to Anak Krakatau and compute the average horizontal wind velocities as a function of height above sea level from 20:00 to 23:00 (Additional file 1).

Detection using fixed-source time-frequency scanning
We develop a simple and fast detection method by reversing the application of existing seismic source localization methods using the spatial distribution of amplitudes at the expected arrival times across a seismic network (Kumagai et al. 2009(Kumagai et al. , 2010)), while also being inspired by similar methods using seismogram phases (Kao andShan 2004, 2007).These methods measure the likelihood of a trial location as the seismic source, eventually searching for an optimum source location in the spatial domain, given seismic signals at a specified time window and frequency band.Our idea is that instead of searching for an optimum location in the spatial domain, we scan the seismograms in the time and frequency domains and measure the likelihood of a seismic event occurring at a fixed, potential location.At an assumed seismic source location with an origin time of t and a frequency of f , we may measure the source likelihood by simply averaging the absolute seismic amplitudes at the expected arrival times at N sta- tions (Kao and Shan 2004): in which u i , τ i , and t + τ i are the seismogram, expected travel time, and expected arrival time at station i , respec- tively.Larger B values represent higher likelihoods that the seismic signals recorded at the stations originated from the assumed source.Here, we account for attenuation by applying simple seismogram normalization (e.g., Kao and Shan 2004).
However, the amplitudes at the expected arrival times may not be accurate due to medium heterogeneities causing seismic energy to be delayed or spread in time, (1) and local amplification may cause deviations in the expected amplitudes.For the detection method to be less sensitive to these uncertainties, we first calculate the seismogram envelope ε i by taking the absolute Hilbert trans- form.Instead of using only amplitude at t + τ i , we use the average of envelope amplitudes in a time window with a width of L and centered at t + τ i (e.g., Kao and Shan 2007): where M is the number of data samples inside the time window.Then, we account for attenuation by assuming the attenuation of surface waves through isotropic radiation in a homogeneous medium (e.g., Battaglia and Aki 2003): where A 0 is the source amplitude, A i is the theoreti- cal amplitude at station i , and 1/ √ r i is the theoretical attenuation from geometrical spreading with epicentral distance r i from the source to station i .The exponential term is the theoretical anelastic attenuation with seismic quality factor Q and phase velocity α .Prior to the scan- ning, we apply an initial correction to each station by adopting the local amplitude correction from Permana and Aoyama (2023): which is analogous to the site amplification factors commonly computed using coda waves of local and regional tectonic earthquakes (Kato et al. 1995;Phillips and Aki 1986).However, s i may not only correct for site ampli- fication near the ground surface where our stations are located.Previous studies have shown that the assumption of isotropic radiation is typically valid at high frequencies of 2-5 Hz or higher (e.g., Maeda and Obara 2009;Takemura et al. 2009), where the seismic radiation pattern is distorted due to scattering from local-scale medium heterogeneities.At lower frequencies, the effect of the radiation pattern may become more pronounced.Haney (2010) measured the vertical amplitudes of surface waves from volcanic tremors below 0.5 Hz over a local seismic network in Alaska and observed the radiation pattern from a particular source mechanism instead of isotropic radiation.Applying the local amplitude correction as: (2) means that theoretically, the corrected amplitude E i should be equal to A i .Therefore, the s i removes the com- bined effect of site amplification and radiation pattern.Due to the uncertainties from medium heterogeneities, and depending on how we estimate the unknown parameter A 0 , E i becomes only an approximation of A i .To increase the confidence of the detection, we suggest estimating A 0 and s i using seismograms of target events that actually originated at the assumed source.
For the scanning over various t and f , we redefine the calculation of B by further correcting E i for the expected attenuation of surface waves as: Expecting the observed amplitude attenuation to follow the theoretical one, which heavily depends on the distance between the source and the stations, provide more constraints to the likelihood than simply normalizing the seismograms.Moreover, B is now equivalent to the estimation of seismic source amplitude from Ogiso et al. (2015).Therefore, for a seismic event that originated at the assumed source, B is physically related to its source amplitude.
Another likelihood measure is introduced by examining if the distribution of source-normalized amplitudes E i /B as a function of r i agrees with the expected attenu- ation as: which is similar to the amplitude residual function (Battaglia and Aki 2003;Ogiso et al. 2015), except that in our approach, we allow C to have varying signs to help us dis- tinguish the signatures from different seismic events in more detail.A C value closer to zero represents a higher likelihood of a seismic event to originate at the assumed source, while a higher (positive) or lower (negative) value may represent the amplitude characteristics of different types of seismic events, as will be investigated later.
Our formulations mean that we may scan the seismograms in time and frequency to search for the signature of possible signals coming from a known source with the attenuation characteristics of surface waves.In principle, only a short duration of seismic data is required, from (5 , t + τ i − L/2 at the nearest station from the assumed source to t + τ i + L/2 at the farthest station, where we define τ i = r i /α .Note that the frequency dependence of α , Q , and s i needs to be determined beforehand. The signature of seismic signals will be characterized and identified by evaluating E i , B , and C at selected t and f values.Scanning in time and frequency over a seismic network is similar to the seismic covariance matrix analysis of Seydoux et al. (2016), which searches for coherent signals while the true source remains unknown.Scanning in time to get information from seismic signals at expected arrival times is also similar to the vespa process in array analysis investigating distant earthquakes (Davies et al. 1971).Our detection method is computationally fast, with the prospect of providing a detailed time-frequency scan using near-real-time data, while also focusing on a potential source location.

Parameter tuning for collapse event detection
The surface wave phase velocity α , seismic quality fac- tor Q , and local amplitude corrections s i are determined by examining the seismic event associated with the flank collapse.For this purpose, we use the seismic data from 20:50 to 21:05 on 22 December 2018 (Fig. 1b).We examine the frequency dependence of α by computing the cross-correlations of normalized seismograms at all station pairs and various frequency bands, then observing the moveout in the delay time of cross-correlation peaks as a function of the residual epicentral distance of each station pair from Anak Krakatau.We find that these moveouts can be explained with velocities of 3-3.5 km/s, especially at 0.02-0.5 Hz (Fig. 2a).Getting stable crosscorrelation peaks at > 0.5 Hz and residual distances of > 200 km is increasingly difficult, indicating that seismic signals at higher frequencies are less coherent across stations, particularly at distant stations.The source mechanism of the mass movement during the flank collapse event has been studied using seismograms at lower frequencies, while higher frequencies of > 0.7 Hz are associated with the eruptive activity (Walter et al. 2019).Considering that the variations are relatively small, we simplify the velocity value and define α = 3 km/s for the detection at 0.02-1 Hz.
A single s i for each station and a single Q are typically used within a fixed frequency band in local (e.g., Morioka et al. 2017;Walsh et al. 2017) or regional seismic networks (Cook et al. 2021;Maeda and Obara 2009).For our detection, we need to examine the frequency dependence of Q and s i .However, it is difficult to assume reasonable Q due to the lack of dedicated studies (e.g., Sharma and Mitra 2018) in the area of our seismic network.We adopt and modify the method of Permana and Aoyama (2023) to estimate Q and s i according to the theoretical attenu- ation of known seismic events.We set t to 20:56:00 on 22 December 2018 to ensure that seismic amplitudes at the expected arrival times are associated with the peak of flank collapse signals.We measure E i in a 10-s window with f of every 0.01 Hz using a 0.02-Hz window.We first determine Q by reducing the differences between the observed amplitudes E i and the theoretical ones A i by minimizing the following root-mean-squared (rms) error: searching over various trials of Q and source ampli- tudes A 0 until finding a minimum R .Epicentral dis- tances r i are computed from Anak Krakatau.We find that the optimum Q that produces the minimum R increases with frequency, which tends to be more stable at < 0.2 Hz and greatly vary at > 0.5 Hz, probably due to the low-coherence surface waves at higher frequencies (Fig. 2b).The more stable Q at lower frequencies up to around 0.1 Hz has been noted by previous studies (e.g., Sato and Fehler 1998).For the detection method, we define an empirical Q estimation by manually fitting a power-law curve (Fig. 2b): We can now estimate the local amplitude correction for each station by repeating the minimization of R (Eq.8) with the above empirical Q estimation, obtain- ing the optimum A 0 , and computing s i .The obtained s i vary significantly with frequency and thus cannot be simplified into a single value.Therefore, we smooth these frequency-dependent s i using a 5-sample moving aver- age (Fig. 2c and Additional file 2).At > 0.3 Hz, the local amplification correction s i may not correctly reflect the amplitude differences between the flank collapse event and the theoretical attenuation, because the collapse event is more dominant at lower frequencies.Particularly ( 9) at distant stations, local noises may affect the resulting s i , such as at MNAI station, where very high s i values at > 0.3 Hz (Fig. 2c) are obtained due to the high level of local noise dominating the flank collapse signals.
In order to see if s i can properly correct for the ampli- tude variations other than the expected attenuation, we compute the locally uncorrected amplitudes E i and the corrected ones E i at 0.02-0.04Hz, 0.08-0.1 Hz, and 0.2- 0.3 Hz, with f , Q , and s i averaged within each frequency band.We normalize those amplitudes with B (equivalent to the source-normalized amplitudes) and compare their spatial distribution with the expected attenuation of surface waves exp −πfr i /αQ(f ) / √ r i assuming an isotropic source at Anak Krakatau (Fig. 3).Although we do not investigate the source mechanism and radiation pattern in this study, we find that E i are more similar to the theo- retical amplitudes than E i , showing that s i is able to com- pensate for the presumed effects of the local site effects and radiation pattern.

Detection application and evaluation
The 2018 flank collapse and tsunami episode With Anak Krakatau as our fixed source, we apply the detection method to scan the seismogram with f of every 0.01 Hz using a 0.02-Hz window and t of every 0.5 s using a 10-s window.The length of the time window has the same effect as smoothing, where longer windows cause E i to be less sensitive to amplitude changes in time, thus reducing the temporal resolution of the scanning.We find that this smoothing effect no longer changes significantly using time windows of 10 s or shorter.We perform the scanning for two sets of stations: those with an epicentral distance of < 400 km (all 13 stations) and < 200 km (4 stations) from the volcano, to examine the effect of the number and distance of the seismic stations.
Events that are visually distinct in the seismograms (Fig. 4a) correspond to high B values at the frequen- cies where signal energy is dominant (Fig. 4b and c).In our detection method, a non-monochromatic seismic event that originated at the assumed source with an α of around 3 km/s should show higher B with onsets that are aligned in time over the frequencies where its seismic energy is distributed.This characteristic is used as the first indication of an event coming from Anak Krakatau, such as the flank collapse event showing rapidly increasing B with sharp onsets at around 20:55.However, the signals interpreted as a teleseismic event at 21:35 also show time-aligned high B .On the other hand, the high B during the following low-frequency signals around 21:55 does not appear to be aligned in time.Therefore, using only B is not enough to distinguish signals coming from Anak Krakatau from the other signals.Note that for those other signals, B does not correctly express the source amplitude, although it still attempts to do so by correcting the amplitudes according to the surface wave attenuation.
Next, we examine if C values can distinguish these sig- nals from each other.We exclude C from the first and last 700 time windows to remove the artifacts due to the effect of tapering and limited data at both ends of the seismograms.The flank collapse event is associated with C of around zero and higher (positive), while the impulsive teleseismic and the following low-frequency signals correspond to negative C values (Fig. 4b and c).
To better understand what C means, we examine three ranges of − 0.3 ± 0.05, 0.0 ± 0.05, and 0.3 ± 0.05 representing negative, zero, and positive C , respectively.From the scanning result using stations at < 400 km epicentral distance (Fig. 4b), we randomly select 1000 time-frequency points within each range, collect the E i , and calculate the source-normalized amplitudes averaged at each station.We also compute the corresponding expected surface wave attenuation curves exp −πfr i /αQ(f ) / √ r i for epi- central distances of 0-400 km and average them over all points.For C ≈ − 0.3, E i tends to increase with epicentral distance, which is the opposite of the expected attenuation behavior (Fig. 5a).This explains the negative C asso- ciated with the teleseismic signals, where the earliest and largest amplitudes arrived at the farthest stations from Anak Krakatau in the direction of incoming signals.The E i values for C ≈ 0 tend to closely follow the theoretical attenuation curves, which is the expected behavior for signals originated at Anak Krakatau (Fig. 5b).However, during the flank collapse event, positive C values are also often obtained.Examination at C ≈ 0.3 (Fig. 5c) reveals higher E i at short epicentral distances and lower E i at longer distances, but nevertheless remains roughly consistent with the expected attenuation from Anak Krakatau.The examination of C is able to further distinguish the flank collapse event, from which the analysis parameters are derived, from the other detections of high B values.
The background values of B are small and persistent throughout our data duration and are associated with seismic noise.The characteristics of noise are particularly clear before the flank collapse, during a period of relative seismic quiescence associated with conduit sealing (e.g., Roman et al. 2016).We observe relatively higher B at 0.14-0.3Hz associated with the secondary microseism from the interaction of wind-driven sea waves and sea waves-seafloor interaction around the coast (Friedrich et al. 1998;Hasselmann 1963).B at higher frequencies of > 0.2 Hz are generally larger than those at < 0.1 Hz, probably due to noise sources more local to the stations.These noises correspond with negative C , suggesting a different explanation than that of the teleseismic signals.
Since noises are less correlated between stations and are more local, their amplitude levels are more similar to each other and do not necessarily follow the expected attenuation from a particular source.Persistent signals started at the onset of the flank collapse at > 0.35 Hz may be related to the eruptive activity triggered by the partial collapse of the volcanic edifice (e.g., Walter et al. 2019) and will be examined later.
The scanning result using four stations at epicentral distances of < 200 km (Fig. 4c) shows similar B char- acteristics to those using all stations but with generally lower values.Nevertheless, B associated with the flank collapse remain the largest.Unsurprisingly, the C values from using fewer stations are significantly more difficult to interpret.This is because a smaller number of E i increases the uncertainty when com- paring them with the theoretical attenuation.Adding more seismic stations, preferably at varying epicentral distances, may produce a more robust value and reduce the uncertainty of C , but with the trade-off of increasing B for seismic signals unrelated to the assumed source including seismic noise, as noted at higher frequencies.Moreover, stations that provide better azimuthal coverage of the assumed source are preferred for better spatial constraints and signal discrimination.In this study, our network of 13 stations provides sufficient coverage in azimuth and epicentral distance among the nearest available stations from Anak Krakatau.

Seismic events not originated from Anak Krakatau
We select four tectonic earthquakes (EQ1-EQ4, Fig. 1a) to further evaluate the detection characteristics of seismic events that did not originate from Anak Krakatau.We use the origin times and hypocenters from the BMKG earthquake catalog shown in Table 1.EQ1 and EQ2 occurred in Sunda Strait, less than 20 km from Anak Krakatau with a focal depth of 51 and 128 km, respectively.Thus, we may examine the cases where seismic signals are coming from near Anak Krakatau at various depths.EQ3 occurred south of western Java, 134 km from Anak Krakatau at 46 km depth.EQ4 occurred about 870 km northwest of Anak Krakatau in the west of Sumatra at 13 km depth, well outside our seismic network with a distance of around 530 km from MNAI station.Therefore, with respect to our seismic network and epicentral distance from Anak Krakatau, these earthquakes represent the cases of local earthquakes (EQ1 and EQ2), earthquakes at intermediate distance (EQ3), and distant earthquakes outside the seismic network (EQ4).The magnitude of these earthquakes ranges from 4.6-6.1,roughly representing the calculated magnitude of the flank collapse event from previous studies (e.g., Walter et al. 2019;Ye et al. 2020).
We scan the seismograms of each earthquake using the same source and parameters as in the case of flank collapse and tsunami episode.The results show that the B associated with the earthquakes are high, but with different C behavior between earthquakes (Fig. 6).For local and intermediate-distance earthquakes (EQ1-EQ3), body waves dominate the waveforms at most stations, particularly for deeper earthquakes, associated with large B at around 0.2 Hz or higher (Fig. 6a-c).Since EQ2 has a smaller magnitude, a deeper source, and therefore a lower SNR, its B values tend to be obscured by the sec- ondary microseism at around 0.2 Hz.Meanwhile, a distant earthquake at a shallower depth (EQ4) is dominated by surface waves at lower frequencies (Fig. 6d).The body waves of EQ1 and EQ2 are associated with high positive C (Fig. 6a-b), while EQ3 does not show a clear distinc- tion with body wave-associated C values of around zero (Fig. 6c).On the contrary, EQ4 is associated with negative C from before the time of high B corresponding to the lower-amplitude body waves, followed by those from the more dominant surface waves (Fig. 6d).
We further examine the distribution of source-normalized amplitudes with epicentral distance from Anak Krakatau within the selected "examination windows" of different C characteristics shown in Fig. 6: body waves from EQ1 and EQ2 ( C > 0.3), body waves from EQ3 ( C ≈ 0), and surface waves from EQ4 ( C < − 0.3).Within each examination window, we collect 1000 time-frequency points that satisfy the corresponding C condition.We average the source-normalized amplitudes for each station and compute the corresponding theoretical surface wave attenuation curves, similar to the previous analysis in Fig. 5.The body waves of local earthquakes correspond to higher E i at nearby stations and lower E i for more dis- tant stations (Fig. 7a), which are similar to the positive C from the analysis of flank collapse and tsunami episode (Fig. 5c) but with a higher extent of deviations from the expected attenuation.This may also be consistent with the higher attenuation of body waves ( 1/r i ) than that of surface waves ( 1/ √ r i ).However, we cannot fit a theoreti- cal attenuation of 1/r i to our E i because our amplitude corrections and parameters have been tuned using surface waves.
For the case of C ≈ 0 during EQ3 (Fig. 7b), the E i dis- tribution does not follow the theoretical attenuation like that for signals originated from Anak Krakatau (Fig. 5b).Note that we do not use all stations in our network due to amplitude anomalies for stations located very close to the earthquake epicenter.Nevertheless, we observe that a value of C ≈ 0 may not always represent seismic sig- nals coming from Anak Krakatau, and therefore further examinations of E i distribution are important to avoid false detections.EQ4 shows the largest E i at stations closest to the epicenter, which are the farthest from Anak Krakatau (Fig. 7c).Similar to the teleseismic signals, the farthest stations from both the epicenter and Anak Krakatau recorded the smallest amplitudes.Despite the fact that both the largest and smallest amplitudes are obtained at the farthest stations, negative C will be more likely to be obtained for earthquakes outside the seismic network, because larger amplitudes affect the C calculation more than the smaller ones.The consistency of negative C for both body and surface waves means that we can easily identify a seismic event that occurred outside our seismic network.
The evaluation of the detection method using data from the flank collapse and tsunami as well as tectonic earthquakes reveals different signatures for seismic events at different epicentral distances from the assumed source.In principle, the stations closest to a seismic source will measure the highest seismic amplitudes.The further that source is from the assumed source (i.e., Anak Krakatau), its E i distribution becomes increasingly differ- ent from the expected attenuation from the assumed one and therefore becomes easier to identify as an event that did not originate at the assumed source.
We do not discuss the origin time of the tectonic earthquakes, since origin time will be correctly measured only for an event that originated at Anak Krakatau.We also ignore C ≈ 0 or highly positive C related to very low B associated with seismic noises.Such C values do not dis- appear when changing the time window length for scanning and may be related to various uncorrelated noises local to the stations.

Seismic signatures during the 2018 flank collapse and tsunami
With the knowledge of the detection behavior for different seismic events, we can now investigate the seismic signatures during the 2018 flank collapse and tsunami episode in more detail.In the following discussions, we examine the scanning results using seismic stations at < 400 km from Anak Krakatau (Fig. 4b).We inspect each seismic signature of interest in an examination window of a specified time and frequency range, within which we collect all available E i , B , and C values.Finally, we com- pute the average C , compare the distribution of averaged source-normalized amplitudes with the epicentral distance from Anak Krakatau with that of the corresponding averaged surface wave attenuation curves (following Figs. 5 and 7), and map the averaged source-normalized amplitudes at the respective station locations.We also average B over frequencies of interest to help determine the timing of seismic signatures.

The flank collapse event
The flank collapse event in particular has been comprehensively discussed by previous studies (e.g., Perttu et al. 2020;Walter et al. 2019;Ye et al. 2020).Its body wave content is noticeably low, as evident from the comparison of our results (Fig. 4) with those of tectonic earthquakes (Fig. 6).B associated with body waves during EQ1-EQ3 are dominant at > 0.3 Hz, while the flank collapse event is dominated by surface waves at 0.05-0.3Hz, as is similar with EQ4.Perttu et al. (2020) determined that the flank collapse event is dominated by Love waves based on the examination of three-component seismograms.
We focus our analysis on the frequency band of 0.02-0.3Hz, where a clear, rapid increase of B relative to the prior noise-related level is observed.We estimate the origin time of this event by averaging B values at 0.02- 0.3 Hz and marking the approximate start time of the rapid B increase, which is around 20:55:30 (Fig. 8a).Wal- ter et al. (2019) determined the origin time of 20:55:49 from seismograms at 0.1-4 Hz and inverted the signals at 0.01-0.03Hz.They determined a magnitude of M w 5.3 with a non-double-couple mechanism indicating tensile opening and shear rupture with about a 61° dip to the southwest, which is consistent with the visible loss of the volcanic edifice on the southwest side.The likely collapse mechanism that can explain the observed tsunami height is thought to be a single flank failure (Grilli et al. 2019;Paris et al. 2020;Priyanto et al. 2021).
The flank collapse event is clearly detected in our results, associated with high B and C ≈ 0 or positive (Fig. 8a), which becomes our benchmark for detecting similar events.Since this event has been well studied, we examine our scanning results to find new details surrounding it.We note that after the collapse event around 20:55, there are two smaller B peaks that are aligned in time around 21:08:30 and 21:15:15 (Fig. 8a) with the same B level as those associated with secondary microseism and therefore can only be distinguished at < 0.14 Hz.We refer to the flank collapse event and these two later peaks as events F1, F2, and F3, respectively.Seismogram examinations show that F2 and F3 can be distinguished at 0.05-0.35Hz but are difficult to recognize in wider frequency bands (Fig. 8b and Additional file 3).Due to the visual similarity in the B signature, we further compare their distributions of source-normalized amplitudes.
We define examination windows of 60 s, 40 s, and 30 s in time and 0.3 Hz, 0.08 Hz, and 0.05 Hz in frequency where F1-F3 are clearly observed, respectively (Fig. 8a).Note that the width of the examination window is smaller in time and frequency for later events, as they are more difficult to recognize in the scanning results.The averaged source-normalized amplitudes for F1-F3 show the same signature of seismic events coming from Anak Krakatau, despite the deviations from the theoretical attenuation being understandably larger for F2 and F3 (Fig. 8c).The averaged C for F2 and F3 are close to zero and tend to be negative due to being closer to noise amplitudes, but the spatial distributions of the averaged source-normalized amplitudes indicate that they originated in the vicinity of Anak Krakatau (Fig. 8d).
F2 and F3 have not been discussed in earlier studies, as discussions were focused on F1 which is associated with the flank failure process that explained the tsunami generation (e.g., Ye et al. 2020).However, a summary of eyewitness accounts from close distances around Anak Krakatau stated that after the main collapse (F1), there were two additional collapses (Perttu et al. 2020).It was not clear whether these later collapses were observed visually or audibly, as Anak Krakatau was also reported to be obscured from view.The interval between the main and second collapses was reported as being longer than that between the second and third collapses.The time scale of these accounts is not clear, but we assume that these reported events occurred within about 30 min of the main collapse, before the tsunami reached the coasts of Java and Sumatra.In our results, the interval between the onset of F1 and F2 is about 13 min, while that between F2 and F3 is about 8 min.Since these timings appear to be consistent with the eyewitness accounts and the characteristics of F1-F3 are similar, especially at low frequencies, we interpret F2 and F3 as two possible much smaller collapse events that occurred after the main collapse.These additional collapse events may hold new insights on the phenomena surrounding the flank collapse and tsunami episode and will be discussed in later sections.(Fig. 9a).While the main collapse signals at lower frequencies (F1) last for about 8 min, those at higher frequencies show larger amplitudes that vary irregularly for about 13 min before gradually decreasing and becoming more stable until the end of our data period.Averaged B values (Fig. 9b) show that after a peak at the onset of the main collapse (referred to as G1), there is a smaller peak just before the seismic activity becomes less intense (referred to as G2), around the time of the possible second collapse (F2).B values from G1 to G2 show a sustained peak at > 0.35 Hz, which later becomes more prominent at 0.45-0.7 Hz.These values are varying and do not appear to be continuous with time, a characteristic that is actually easier to recognize when using seismic stations at < 200 km from Anak Krakatau (Fig. 4c).These B values are associated with C of around zero or positive, suggesting their point of origin at Anak Krakatau.Note that we ignore the very high C that are likely to be an arti- fact around the arrival of the teleseismic signal.
The distributions of averaged source-normalized amplitudes are examined using 60-s examination windows representing G1, G2, and another window (G3) around 21:48 after the teleseismic arrival (Fig. 9b).We confirm that the measured amplitudes during G1-G3 follow the expected attenuation from Anak Krakatau, with an average C of around zero (Fig. 9c).Their spatial pattern also suggested a source at Anak Krakatau (Fig. 9d).
Seismic signals from the flank collapse at higher frequencies have been associated with volcanic eruptions, as noted by Walter et al. (2019) at 0.7-4 Hz.Perttu et al. (2020) analyzed the spectrograms at several stations and found spectral lines at < 2 Hz, which they interpreted as the repetitive occurrences of explosive sources every 30-50 s associated with Surtseyan eruptions.The noncontinuous B and C in our analysis may correspond with such a source process.The signature of eruptive activity in our results appears to extend to > 1 Hz, beyond the range of our investigation.Therefore, we cannot provide a full interpretation of the eruptive activity only from our analysis at 0.35-1 Hz.The more intense activity from G1 to G2 may represent the more violent stage of phreatomagmatic eruptions due to sudden depressurization of the conduit and contact with a large amount of seawater (Hunt et al. 2018;Williams et al. 2019) or the sudden enlargement of the vent (e.g., Eibl et al. 2023).Eventually, a more stable state of interaction between volcanic fluids (magma and gas) and seawater was reached, producing sustained phreatomagmatic activity.

The teleseismic event
Our initial information regarding the teleseismic event is that the impulsive onset at 21:35 was associated with an earthquake from Vanuatu (Perttu et al. 2020).We search the earthquake catalog of the United States Geological Survey (USGS, https:// earth quake.usgs.gov/ earth quakes/ search/) on 22 December 2018 and confirm that an M w 6.0 earthquake had occurred near Vanuatu at 21:25:01 WIB, at a 6769 km distance and 101.9° azimuth from Anak Krakatau with a focal depth of 42 km (Additional file 4).The signal arrivals at 21:35 are first observed at our easternmost CMJI station and propagated westward (Fig. 10a).However, upon examination of the seismograms at longer duration, we note three distinct signals that show the same characteristics of coming from the east direction and are associated with high B at different frequencies: the 21:35 arrivals at 0.13-0.7 Hz, the signals arriving 7 min later around 21:42 at 0.08-0.25 Hz, and the following low-frequency signals arriving around 21:55 at 0.03-0.07Hz (Fig. 10b).We refer to these three signals as sub-events T1, T2, and T3, respectively; all are associated with negative C .Note that T3 signals can be clearly observed from the raw seismograms, but their emergent onset hinders an accurate reading of their arrival time, while T2 signals are not clearly observed at many stations, especially in wider frequency bands (Fig. 1b).The arrival times of T1 and T2 are inconsistent with surface wave propagation from the earthquake epicenter, but those of the large amplitudes during T3 can be explained with wave velocities of 2.8-3.8km/s, which is approximately similar to our α .By treating our stations as a seismic array due to the very large distance from the earthquake epicenter, we estimate a very high apparent velocity of 15 km/s for T1 and 12 km/s for T2, indicating that T1 and T2 are likely to be the body waves of the teleseismic earthquake, while T3 is the following surface wave phases.
The distribution of the averaged source-normalized amplitudes for 60 s examination windows during T1-T3 (Fig. 10b) shows that the largest amplitudes are always observed at the stations located east of Anak Krakatau (Fig. 11a and b).Note that E i may not always be the larg- est at the easternmost station due to the effect of s i cor- rections.We only examine T3 at 0.03-0.05Hz when B is the most aligned in time (Fig. 10b).The average C values are < − 0.3 for T1 and T2, and < − 0.17 for T3, indicating the signature of a seismic event located outside of our seismic network.
We further seek confirmation that T1-T3 correspond to the same earthquake by calculating the theoretical arrival times of various seismic phases using an implementation of the TauP toolkit (Crotwell et al. 1999) with the iasp91 earth model (Kennett and Engdahl 1991) and the origin time and hypocenter from the USGS catalog.In Fig. 11c, the seismogram at CMJI station is filtered at several frequency bands, showing signal arrivals that coincide with the theoretical arrival times of several P and S-wave phases, with T1 and T2 mainly corresponding with direct P and S waves, respectively.Similar comparisons at all stations are provided in Additional file 4. The larger amplitudes of T3 become more delayed with increasing frequencies, demonstrating the dispersive characteristics of surface waves where velocity decreases with increasing frequency (e.g., Levshin et al. 1989).This surface wave dispersion of teleseismic signals explains why the associated B values do not align in time (Fig. 10b).
Therefore, our findings are not in agreement with those of Perttu et al. (2020), who associated the low-frequency signals (T3) with the tsunami arriving at the coast.Instead, we associate these signals with the surface wave phases of the same teleseismic event that also generated the signals arriving at 21:35 (T1).

Detection of weak teleseismic signals
The detection of the teleseismic event from Vanuatu has the advantage of the visually distinct signals that are clearly observed at many seismic stations.In this section, we show that our detection method is also capable of capturing the signature of another teleseismic event that is difficult to distinguish from the seismograms.Our first clue is a weak B increase at < 0.1 Hz between the time of F2 and T1 in the previous analyses.Filtering the seismogram at 0.04-0.09Hz reveals the amplitudes larger than F2 and F3 that appear to be first arriving at the PMBI station (Fig. 12a).The corresponding B val- ues are lower than those of the secondary microseism, but they show the dispersive characteristics of surface waves similar to T3 (Fig. 12b).We search for possible body wave detections prior to these surface wave signals, but the seismograms were dominated by the signals from the collapse events (F1-F3).Nevertheless, around 20:40, we find a particular weak B increase accompanied by negative C at 0.06-0.1 Hz, which corresponds to sig- that also appear to be arriving first at PMBI with a very high apparent velocity.Furthermore, the averaged B at 0.04-0.09Hz shows a small peak around 20:51 which may be related to another body wave signal.We refer to these weak B peaks at 20:40 and 20:51 and the surface waves between F2 and T1 as sub-events K1, K2, and K3, respectively.We consult the USGS earthquake catalog for a teleseismic earthquake that can explain the first arrivals at PMBI station, preferably in the north direction, and find the most potential M w 5.6 earthquake occurring east of Kamchatka Peninsula at 20:29:46 WIB, with an 8691 km epicentral distance, a 30.2°azimuth from Anak Krakatau, and a 10 km focal depth (Additional file 5).This earthquake is able to explain the arrival times of K3 using velocities of 2.2-3.2km/s (Fig. 12a).
The 60 s, 40 s, and 60 s examination windows for K1-K3, respectively (Fig. 12b), show average C values of < − 0.2, but a clear indication of teleseismic origin from the source-normalized amplitudes is only observed during K1 (Fig. 13a).The amplitudes during K2 and K3 are more similar between stations, revealing a signature similar to seismic noise, probably due to low SNR.The spatial amplitude distribution for K1 shows the largest amplitude of K1 at PMBI and JCJI, indicating signals arriving from the north-northeast direction, while K2 and K3 do not show such an indication (Fig. 13b).Through comparison with the theoretical arrival times of body wave phases, we confirm that K1 is consistent with the arrival of P waves (Fig. 13c and Additional file 5).Possible signals related to P-wave arrivals may also be observed at > 0.5 Hz, but they are difficult to recognize using B and C values (Fig. 9a and b).K2 appears to be consistent with the arrival of S waves, but it is difficult to discriminate the associated signals from the seismograms.Moreover, it is generally very difficult to recognize this earthquake from the seismograms due to the dominating signals of the collapse events and the teleseismic event from Vanuatu.Nonetheless, our scanning results are able to discriminate this earthquake by examining the signal amplitudes at various frequencies over many seismic stations.

Timeline during the flank collapse and tsunami episode
We summarize our findings and construct a simplified timeline to understand the sequence of events during the 2018 Anak Krakatau flank collapse and tsunami episode.We only focus on the time period of our seismic data, discussing various phenomena during the episode and combining them with results from previous studies.We use the timing of the seismic events found in this study, the origin time of teleseismic earthquakes from the USGS catalog, and tsunami arrival times from Grilli et al. (2019) and Paris et al. (2020).To determine the details of related surface phenomena, we use insights from Doppler radar observation.
Because of the 135 km distance from the radar and the available elevation angles, we can only observe from a minimum height of about 4.6 km over Anak Krakatau (Fig. 14a).We focus on the area around the volcano and compute the maximum reflectivity factors in vertical and latitude directions to capture the time evolution of the eruption plume horizontally and vertically, respectively (Fig. 14b).The reflectivity factor has been shown to be proportional to the particle size and the amount of tephra in volcanic plumes and has been used in quantitative volcanic ash identification (Marzano et al. 2006).However, due to the lack of similar studies, we will interpret our radar observations qualitatively.
Around the time of the devastating episode, a persistent, localized strong contrast of reflectivity factors was observed above Anak Krakatau and was stable over time (Fig. 14b).We confirm that this is not a meteorological feature (e.g., cloud or precipitation) by examining the radar data over a longer time period, during which meteorological clouds can be seen passing over the volcano over time with evolving shapes (Additional file 6).Such meteorological characteristics were not observed around Anak Krakatau during the time period of our seismic data, and therefore we interpret these localized reflectivity factors as the eruption plume.During this episode, the wind at 5-20 km heights was stable and blowing in westsouthwest directions (Additional file 1).
Our timeline of events is shown in Fig. 15.The eruptive activity on 22 December 2018 has been active since several hours prior to the flank collapse, with lava flow becoming increasingly intense, followed by more frequent ash emission and the formation of volcanic plumes (Paris et al. 2020;Perttu et al. 2020).The radar captured the start of a stable plume at about 5 km height around 19:00, which shortly disappeared around 20:00 but soon reappeared.Within about 25 min before the flank collapse, despite the fact that the seismicity at regional stations below 1 Hz was relatively quiet, the eruption plume rose below 7 km height (Fig. 14b).Within 10 min before the flank collapse, we detected weak arrivals of the body waves from a M w 5.6 earthquake in the east of the Kamchatka Peninsula.Walter et al. (2019) found a high-frequency signal at > 1 Hz about 115 s before the collapse, which was interpreted as a possible precursor of the flank failure.However, this signal was only clearly observed at the nearest station CGJI, and its source was not well understood.Perttu et al. (2020) suggested that this signal was generated at the failure plane, but with some inconsistencies when compared with infrasound signals.Considering that the timing of this signal coincides with the arrivals of the teleseismic S waves, we suggest another possibility that this high-frequency signal is one of the S-wave phase arrivals (Fig. 13c).
The major flank collapse at around 20:55 generated various observed phenomena: the collapse of the southwest sector of the volcano, the following generation of the tsunami, the sudden increase in eruptive activity due to rapid changes in the conduit and conditions, and subsequently, the rapid increase in the eruption plume height.The seismic signals at < 0.35 Hz have been associated with the collapsing mass movement, while tremorlike signals at higher frequencies have been attributed to the intense eruptive activity during the first 13 min after the collapse (Fig. 9a).Following the collapse, the eruption plume rapidly extended vertically to about 16 km height within less than 10 min, suggested to be dominant in ash content (Gouhier and Paris 2019).The increasingly stronger reflectivity factor at the lower part of the plume is a typical feature found in radar observation of volcanic plumes, representing gravitational particle sorting where During 21:04-21:20, radar observation reveals that the plume height rapidly decreased to about 9 km, despite satellite-based observations capturing a continuous plume probably related to water vapor, volcanic gases, and ash (e.g., Gouhier and Paris 2019).Its horizontal extent became broader and elongated in the southwest direction, consistent with the wind direction and possibly affected by the direction of the lateral decompression during the flank collapse (Gouhier and Paris 2019).However, it is likely that the decrease in plume height was caused by the internal process and/or the eruptive process.The reflectivity factors were also gradually decreasing with time, indicating the decreasing intensity of eruptive outputs.We interpret that the resulting decrease in the buoyancy force leads to the gravitational collapse of the eruption plume, particularly for larger volcanic particles.The extent of this collapse at the surface is not known, but it may help in further obscuring Anak Krakatau from visual observation, as reported around 21:16 (Paris et al. 2020;Perttu et al. 2020), while the wind further distributes the airborne volcanic particles in the southwest direction.It is also not clear if there was a generation of significant pyroclastic flow following the plume collapse, as was proposed to explain the previous 1883 eruption of Krakatau (Deplus et al. 1995).
Around 21:08, we detected signals with a signature similar to the flank collapse but much weaker, which we interpreted as the possible smaller second collapse (Fig. 8).The short time interval between this collapse and the following decrease in eruption activity may indicate a related source process.Following the main collapse which left out material instabilities (e.g., Hunt et al. 2018), the second collapse may have altered the vent condition further, making it more stable and eventually producing a more sustained stage of phreatomagmatic eruption, as indicated by the stable seismic signature at 0.45-0.7 Hz (Fig. 9b).During the horizontal broadening of the collapsed plume around 21:15, we suggest that another possible small flank collapse has occurred.The plume collapse may increase the material instabilities in the already unstable remaining volcanic edifice, triggering additional collapse events.Further triggers for these smaller flank collapses may have also come from the eruptive activity (Hunt et al. 2018).These latter collapse events may be consistent with those mentioned in the eyewitness accounts (Perttu et al. 2020).Their intervals of occurrences, however, are not consistent with Williams et al. ( 2019), who suggested multiple flank failures followed by phreatomagmatic eruptions in rapid succession as the mechanism that removed the southwest flank of Anak Krakatau.Finally, radar observation from 21:36 until the end of our seismic data period suggested that the plume had fully collapsed, with plume height returning to 5-7 km or the same height prior to the main flank collapse.Later, from 22:00 afterwards, the observed plume gradually intensified, and the phreatomagmatic activity remained high until the following days (e.g., Gouhier and Paris 2019;Prata et al. 2020), as shown in Additional file 7. The resulting continuous volcanic plume has been suggested to have produced a freezing effect and electrification in the upper troposphere, resulting in a volcanic thunderstorm (Prata et al. 2020).
From around 21:24, the tsunami was estimated to have arrived on the western coast of Java, followed by arrivals on the southern coast of Sumatra from 21:34 (Paris et al. 2020).We do not detect particular seismic signals or signatures that may be related to the tsunami.Weak signals detected at 0.04-0.09Hz during the propagation and arrival of the tsunami are confirmed to be the surface waves of the M w 5.6 teleseismic earthquake (Fig. 12).Then, the seismograms are dominated by a M w 6.0 teleseismic earthquake from Vanuatu.The corresponding P, S, and surface waves are clearly detected at different frequency bands (Fig. 10), where the surface waves were previously interpreted as the tsunami arriving at the coast by Perttu et al. (2020).However, the corresponding signals arrived about 15-25 min after the observed and estimated tsunami arrival times from Grilli et al. (2019) and Paris et al. (2020).The fact that both the teleseismic waves and the tsunami arrived first in Java and later in Sumatra may have caused the misidentification of these signals.For the tsunami to be recorded as ground vibrations, the possible mechanism involves the interaction between ocean waves and those with the seafloor, which are consistent with the generation mechanism of microseism (Hasselmann 1963).Short-duration changes in microseismic activity are possible, for example, due to storm activity (Nishida 2017).However, compared with ocean waves due to storm activity, tsunami waves may not be energetic or sustained enough to excite such microseism.

Applicability of the method
We have shown the potential of our detection method to reveal the details of a major episode of activity by examining the time-frequency seismic signatures measured at a known potential source.As a general procedure, a fixed source must be assumed, and the phase velocity α , the seismic quality factor Q , and the local amplitude cor- rections s i must be computed beforehand, as well as their frequency dependence.In our study, we use the surface wave-dominated flank collapse event to specifically tune these parameters for future detections of similar events.Alternatively, local and regional earthquakes may also be used to obtain the site amplification factors (e.g., Kato et al. 1995) as a component of s i , study the variation of Q (e.g., Sharma and Mitra 2018), and compute seismic velocities or perform tomographic studies (e.g., Afif et al. 2021).The selection of the window lengths and shifts in time and frequency depends on the desired resolution and/or computation time and resources, which can be optimized through trials.
We have compared the effects of different numbers of seismic stations on the scanning results (Fig. 4).Localization methods using alternative approaches (e.g., Kao and Shan 2004;Kumagai et al. 2009) often require a minimum of three stations ( N = 3) to properly constrain the seismic source location.We suggest finding an N that produces reasonably clear discrimination between the noise, expected events, and other events, as shown in Fig. 4b.It is important to have these stations distributed at various epicentral distances from the assumed source to reduce the uncertainty of E i , B , and C with respect to the theoretical attenuation.Sufficient azimuthal coverage is also suggested to better discriminate events that do not originate from the assumed source, especially distant and teleseismic earthquakes.
The proposed detection implies that all observed signals must propagate from a single source, which is not always the case in real-world applications where signals from multiple sources with varying magnitude and location may exist in a time window.In our case, the eruption-related signals after main flank collapse were overlapped and dominated by the teleseismic P waves (T1), resulting in negative C followed by positive-value artifacts (Fig. 9b).Teleseismic S waves may have also interfered with the signals from the main flank collapse (F1), but with no significant effects due to their low SNR.Note that low-SNR signals tend to show a signature similar to that of seismic noise.Signals from multiple sources will likely deviate the E i distribution further from the expected attenuation.Multiple distant and teleseismic events will result in more pronounced negative C values due to their first arrivals at the outermost stations.

Practical application for real-time monitoring
We have demonstrated the application of the detection method by scanning the available long-duration seismograms.However, the outlook of this method is to detect seismic events as real time as possible.Therefore, we propose a procedure for near-real-time detection for use in, for example, volcano monitoring.Considering an established continuous seismic data stream and the minimum required data of up to t + τ i + L/2 at the farthest station, a station at a greater distance from the assumed source will cause a longer waiting time to obtain sufficient data.Longer data may be necessary to avoid the taper effect and filtering artifacts, while those before t are already available from stored seismograms.Using a four-core processor of 3 GHz speed and analysis code written in MATLAB, our computation time for one t and 99 frequency values (within 0.02-1 Hz) is in the order of 1 × 10 -5 s, assuming all data and parameters have been loaded to computer memory.
The resulting E i , B , and C are then fed to an algorithm to detect and identify potential events that originated at the assumed source.We may use B averaged over a fre- quency band of interest for initial detection of rapidly increasing values, which can be improved by introducing a threshold above which the corresponding time is then marked as the possible origin time (Fig. 16a).However, this approach will be more accurate for events with an impulsive onset than those with an emergent onset, in which case delayed detections may be made by identifying significant peaks in the averaged B .B can also be used to obtain information about the magnitude of an event that originated at the assumed source, since it is physically related to the projected amplitude at the assumed source.Kumagai et al. (2013) used the source location method utilizing seismic amplitudes (equivalent to E i ) and found a linear relationship between the logarithm of source amplitudes (equivalent to B ) and the magnitude of various volcanic earthquakes.
We have shown that C may faces difficulties in discrim- inating events, such as in the case of C ≈ 0 which may represent an event originated at the assumed source or an intermediate-distance event (Figs. 5b and 7b).Therefore, we propose to compute the similarity between the shape of the source-normalized amplitude curve, as a function of epicentral distance from the assumed source, and that of the theoretical attenuation.We introduce a simple identification scheme by comparing C and this curve similarity, which referred to as the γ function (Fig. 16b).γ can be developed from some known similarity measures, such as the correlation coefficient or cosine similarity.This scheme also involves the selection of three threshold values: Ŵ 1 , Ŵ 2 , and Ŵ 3 .Ŵ 1 and Ŵ 2 are used to discriminate C values into three categories: negative C , C ≈ 0 , and positive C , while Ŵ 3 is a γ limit that is mainly introduced to avoid ambiguous interpretations of C ≈ 0 .In addition, Ŵ 3 is also useful to distinguish seismic noise from distant and teleseismic events for negative C (Figs. 5a and 7c) and identify body waves from local earthquakes for positive C (Figs. 5c and 7a).An event that originated at the assumed source should satisfy Ŵ 1 < C < Ŵ 2 and γ > Ŵ 3 .
The entire detection and identification procedure should be able to provide information about the significance of an event within several minutes after its occurrence, which is then combined with other information in the decision-making process of an early warning system.

Prospects for future developments
A rapid detection and warning system requires many reliable information and insights as earliest as possible for quick decision-making.To support such urgencies, monitoring methods and data at the highest temporal resolution available are preferred.We show the potential of complementary eruption plume monitoring using the available network of weather radars, with a temporal resolution higher than many satellite-based monitoring approaches.Combining an already established regional seismic and weather radar monitoring network with the existing detection system may improve the detection of major events, such as tsunamigenic eruptions or the flank collapse of a volcanic edifice.
The potential of our detection method is not limited to those demonstrated in this study.The detection parameters may be tuned using simulated seismic events, especially in the sense of hazard mitigation and recreation of historical events.Further applications may employ the horizontal components of the seismograms as well as method extensions to detect body wave-dominated events.In such cases, the effect of radiation pattern may be more pronounced, as shown using horizontal seismograms of S waves from a regional seismic network by Takemura et al. (2009) and surface waves from a local network by Haney (2010).Instead of assuming isotropic radiation, the radiation pattern from an expected source mechanism (e.g., Kumagai et al. 2010) may provide more constraints for a more selective detection.Multi-component seismograms may also add more constraints on the polarization of seismic waves.Nevertheless, these details on the target events need to be investigated beforehand for more accurate and reliable detections.

Conclusions
We developed and evaluated a fast detection method for surface wave-dominated seismic events, motivated by the lack of body wave-based detection during the 22 December 2018 Anak Krakatau flank collapse and tsunami.With Anak Krakatau volcano as a fixed seismic source, we measured the seismic amplitudes at the expected arrival times at seismic stations and compared them with the theoretical attenuation of surface waves.By analyzing seismic signatures over a regional seismic network during the flank collapse and tsunami episode, we were able to distinguish seismic signals originated at Anak Krakatau from other events using examples of known tectonic earthquakes.We succeeded in detecting the flank collapse event, from which the method parameters were derived, and the teleseismic event identified in previous studies.In addition, we found other distinct seismic signatures associated with two possible smaller collapse events, changes in the eruptive activity shortly after the main collapse, and the body and surface waves of two teleseismic earthquakes.Combining these results with previous studies and insights on the eruption plume from weather radar observation that revealed volcanic plume collapse shortly after the main flank collapse, we constructed a more detailed timeline around the time of the collapse and tsunami episode.With the potential to detect seismic signatures for different types of events, we further proposed a procedure for real-time event detection and identification from a known potential seismic source that may assist the monitoring of future volcanorelated surface tsunamigenic events.reviewers for their consideration and reviews that greatly improved this manuscript.TY was supported by JST SPRING, Grant Number JPMJSP2119.ASH was supported by the International Joint Graduate Program in Earth and Environmental Sciences (GP-EES) of Tohoku University and JST SPRING, Grant Number JPMJSP2114.

Fig. 1
Fig. 1 Location and seismic episode of interest in this study.a Map showing Anak Krakatau location, distribution of seismic stations (yellow circles), and weather radar location (green diamond).Stars denote the epicenters of selected tectonic earthquakes EQ1-EQ4.Dashed circles illustrate epicentral distances from Anak Krakatau.Insets show detailed location of Anak Krakatau (top right) and the location of EQ4 relative to the stations (bottom left).b Vertical seismograms during the 2018 flank collapse and tsunami from four stations at various epicentral distances

Fig. 2
Fig. 2 Parameter determination for flank collapse event detection.a Seismogram cross-correlations at four frequency bands.Red lines denote the fitted phase velocity to explain the peak delay times.Red circles highlight the peak delay times at higher frequencies.b Measured and fitted seismic quality factors.c Frequency-dependent local amplitude corrections for each station

Fig. 3
Fig. 3 Comparison of local amplitude corrections s i at different frequency bands.Circle colors represent the logarithm of observed amplitudes without corrections E i (top row) and the corrected ones E i (bottom row).Background colors show the theoretical amplitudes (theoretical attenuation) assuming isotropic radiation of surface waves from Anak Krakatau

Fig. 4
Fig. 4 Time-frequency scanning of signals during the 2018 Anak Krakatau flank collapse and tsunami episode.a Seismogram at CGJI station.The scanning is performed using stations with epicentral distances of b < 400 km and c < 200 km from the volcano

Fig. 5
Fig. 5 Source-normalized amplitudes with epicentral distance for the 2018 Anak Krakatau flank collapse and tsunami episode.We compare the amplitudes with a C ≈ − 0.3, b C ≈ 0, and c C ≈ 0.3, from the scanning using stations at < 400 km from Anak Krakatau, whose names are shown in b.Red lines and areas show the average and ± 1 standard deviation of the theoretical surface wave attenuation curves

Fig. 7
Fig. 7Source-normalized amplitudes with epicentral distance from Anak Krakatau for the earthquakes EQ1-EQ4.We compare the amplitudes within the examination windows in Fig.6with a C > 0.3 (from EQ1 and EQ2), b C ≈ 0 (EQ3), and c C < − 0.3 (EQ4).Amplitudes from stations closest to the earthquake epicenter are labeled.Red lines and areas show the average and ± 1 standard deviation of the theoretical surface wave attenuation curves

Fig. 8
Fig. 8 Signature of flank collapse-related events F1-F3.a Scanning results and the identification of F1-F3 examination windows used for computing (c) and (d).Red-shaded areas denote the time length of the windows.b Filtered seismograms at three stations showing F1-F3 signals.c Averaged source-normalized amplitudes with epicentral distance for the examination windows in a. Colors are scaled between the minimum and maximum amplitudes for each event.Red lines denote the averaged surface wave theoretical attenuation.d Maps of the averaged source-normalized amplitudes at the locations of seismic stations

Fig. 9
Fig. 9 Signature of eruption signals related to the flank collapse event.a Filtered seismograms at CGJI and PMBI stations.b Scanning results and the identification of G1-G3 examination windows.c Averaged source-normalized amplitudes with epicentral distance for G1-G3.d Maps of the averaged source-normalized amplitudes at each station

Fig. 10
Fig. 10 Signature of the M w 6.0 teleseismic event from Vanuatu region.a Filtered seismograms sorted according to epicentral distance.Red and blue-dashed lines denote the arrivals of T1 and T2 signals, respectively.Green-shaded area denotes the theoretical arrival times for velocities of 2.8-3.8km/s computed from the epicenter and origin time from the USGS catalog, which coincide with the arrival of T3 signals.b Scanning results and the identification of T1-T3 examination windows for the analyses in Fig. 11a and b

Fig. 11
Fig. 11 Examination of sub-events T1-T3 identified in Fig. 10b.a Averaged amplitudes with epicentral distance from Anak Krakatau for T1-T3 windows.b Maps of the averaged source-normalized amplitudes at each station.c Filtered seismograms from the CMJI station and the theoretical arrival times of body wave phases

Fig. 12
Fig. 12 Signature of the M w 5.6 teleseismic event from east of the Kamchatka Peninsula.a Filtered seismograms identifying K1 and K3 signals.Red dashed line denotes the arrival of K1 signals.Green-shaded area denotes the theoretical arrival times for velocities of 2.2-3.2km/s computed from the epicenter and origin time from the USGS catalog.b Scanning results and the identification of K1-K3 examination windows for the analyses in Fig. 13a and b

F1Fig. 13
Fig. 13 Examination of sub-events K1-K3 identified in Fig. 12b.a Averaged source-normalized amplitudes with epicentral distance from Anak Krakatau for K1-K3 windows.b Maps of the averaged source-normalized amplitudes at each station.c Filtered seismograms from the PMBI station and the theoretical arrival times of body wave phases

Fig. 14
Fig. 14 Time evolution of radar reflectivity factor over Anak Krakatau.a Illustration of the radar scanning planes at various elevation angles (blue lines).Only elevation angles of > 0.5° are used.Red-shaded rectangle represents the vertical area examined in b. b Vertical (upper panel, facing north) and horizontal (bottom panel) distributions of maximum reflectivity factors at each radar observation.White areas denote no radar data.White lines denote the outline of Anak Krakatau and the surrounding islands Fig. 15 A simplified timeline of the 2018 Anak Krakatau flank collapse and tsunami episode.This timeline combines new results from this study with previous studies.Different colors represent different phenomena.Vertical dotted line denotes the time of the major flank collapse

Fig. 16
Fig. 16 Illustrated procedure for real-time event detection.a Initial detection using averaged B .b Identification scheme to distinguish seismic events using C and γ values, along with the thresholds Ŵ 1 , Ŵ 2 , and Ŵ 3 .Green rectangle denotes the preferable result where Ŵ 1 < C < Ŵ 2 and γ > Ŵ 3 .Colors are illustrative and do not represent actual C values

Table 1
List of selected tectonic earthquakes from the BMKG earthquake catalog