Enigmatic Tsunami Waves Amplified by Repetitive Source Events Near Sofugan Volcano, Japan

On 8 October 2023, mysterious tsunamis with a maximum wave height of 60 cm were observed in Izu Islands and southwestern Japan, although only seismic events with body‐wave magnitudes mb 4–5 have been documented to the west of Sofugan volcano. To investigate the source process, we analyze tsunami waveforms recorded by an array network of ocean bottom pressure gauges. Stacked waveforms of pressure gauge records suggest recurrent arrivals of multiple wave trains. Deconvolution of the stacked waveforms by tsunami waveforms from an earlier event revealed over 10 source events that intermittently generated tsunamis for ∼1.5 hr. The temporal history of this sequence corresponds to the origin times of T‐phases estimated by an ocean bottom seismometer and of the seismic swarm, implying a common origin. Larger events later in the sequence occurred at intervals comparable to the tsunami wave period, causing amplification of later phases of the tsunami waves.

• On 8 October 2023 (UTC), enigmatic tsunamis up to ∼60 cm were recorded along broad Japanese coasts without large earthquakes • Analysis of stacked tsunami waveform data shows >10 events intermittently caused tsunamis for ∼1.5 hr with T-phase and seismic excitation • Larger events later in the sequence occurred with intervals similar to the wave period, amplifying the tsunami waves in the later phase

Supporting Information:
Supporting Information may be found in the online version of this article.et al., 2020;Paris, 2015), or submarine/coastal landslides or mass failures (e.g., the 1998 Papua New Guinea tsunami) (Synolakis et al., 2002;Tappin et al., 1999).
This study investigates the source process of the enigmatic tsunami waves using tsunami waveform data of spatially dense ocean bottom pressure (OBP) gauges off the southwestern coast of Japan.We first apply a waveform stacking technique to the multiple tsunami waveforms, and then estimate the temporal history of the tsunami generation based on the analysis of the stacked data.Consequently, we propose their peculiar tsunami origin by repetitive source events that took place intermittently for ∼1.5 hr.

Data
The tsunami waves were recorded by the OBP gauges of Dense Oceanfloor Network system for Earthquakes and Tsunamis (DONET) off the southwestern coast of Japan (Figure 1b) (National Research Institute for Earth Science and Disaster Resilience, 2019).In Figure 2a, we show the OBP data after removing the tidal component, demonstrating repetitive strong high-frequency (>1 Hz) signals approximately from 19:00 to 21:30 (see Figure S1 in Supporting Information S1).These high-frequency signals are T-phases, seismic waves converted from oceanic acoustic waves (Okal, 2008), based on their arrival times explained by a typical T-phase speed of 1.5 km/s from the origin times and locations of the seismic events.Following the T-phase signals, tsunami waves with longer periods were recorded, as shown evidently in the band-pass (0.00125-0.02Hz/50-800 s) filtered records (Figure 2b).Smaller oscillations start around 20:40, leading to the largest amplitudes of ∼20 mm after ∼22:00, exhibiting long-lasting tsunamis for hours with late arrivals of the largest waves.
To capture the features, a tsunami waveform stacking technique is applied to two sets of >10 OBP records of nearby stations: (a) 16 records of DONET1 (KMA-KME), and (b) 13 records of DONET2 (MRC-MRE) (Table S2 in Supporting Information S1).By assuming a point source near the source region of the seismic event swarm on 8 October 2023 (140.026°E,29.787°N), the tsunami travel times to the stations are computed by a shallow-water-wave tsunami model of the Geoware TTT Software (Figure S2 in Supporting Information S1).
The band-pass filtered OBP waveform at each station is shifted as the arrival times are aligned with that at the earliest-arrival station (KMC21 for DONET1, and MRE20 for DONET2).We then stack the time-shifted waveforms and take the amplitude average and the standard deviation.
We show the waveform stacking results with DONET1 in Main Text (Figures 2c-2e) and contain those with DONET2 in Figure S3 in Supporting Information S1.As shown in Figures 2c and 2d, all the DONET1 waveforms here show similar shapes; thereby, the waveform stacking yields clear tsunami waveforms only with small standard deviations.The tsunami waves initiate at ∼120 min (∼20:40) and reach the maximum amplitude at ∼220 min (∼22:10).For further investigation, we apply the wavelet analysis to the stacked waveform (Figure 2e).The obtained scalogram shows that the tsunami signals are composed of multiple bends of dispersive amplitude peaks with early arrival of lower-frequency amplitude followed by higher-frequency amplitude (red arrows and a bracket in Figure 2e).The same feature can be seen in the stacke waveform of DONET2 (Figure S3 in Supporting Information S1).This feature with multiple bands is quite different from a tsunami event originating from a single volcanic earthquake at Sumisu Caldera, which exhibits only a single band (Figure S4 in Supporting Information S1; see the caption for details).Therefore, we speculate that multiple tsunami wave trains, each with a strongly dispersive character, recurrently arrived at the OBP gauges.

Estimation of Tsunami Source Time Function
Using the stacked tsunami waveforms obtained above, we investigate the temporal history of the tsunami generation process using the iterative deconvolution method, which has been widely applied in earthquake source studies.We separately analyze the two stacked waveforms of DONET1 and DONET2 in the same way, but the methodology is explained below with the waveform of DONET1.
We hypothesize that multiple source events took place at similar locations but at different timings, and that each single event produced tsunami waveforms with the same shapes and different amplitudes.Under this hypothesis, the stacked OBP tsunami waveform is the convolution of the temporal history of multiple source events, or the tsunami source time function (STF), and the tsunami waveform produced by a single event, or the Green's 10.1029/2023GL106949 3 of 11 function.Denoting the stacked tsunami waveform as d(t) and the Green's function as w(t), the convolution can be expressed, as follows: where m i and t i represent the source amplitude and timing of the ith iteratively determined source event, respectively, as explained in detail later.
The Green's function is extracted from the stacked tsunami waveform data; in other words, we prepare an empirical Green's function (e.g., Y. Wang et al., 2022).We first confirm that the first seismic event Se01 did not produce measurable tsunami waves, but that the theoretical arrival times of the tsunamis caused by Se02 and Se03, agree well with the timings when the tsunami signal initiates (∼20:40) and when the amplitude increases (∼21:00), respectively (arrows in Figure 3a).We assume that the signal between the two tsunami arrival times represent the tsunami waveform due to Se02, and construct the Green's function by adding an initial zero-amplitude data for the length of the theoretical tsunami travel time to the stacked waveform in the   S2 in Supporting Information S1).The time is aligned with the tsunami arrival at KMC21 (Figure 1b).(d) The stacked waveform (average amplitude; thick black line) ± standard deviation (shaded by gray).Wavelet analysis by the continuous wavelet transformation in MATLAB (Lilly, 2017).Red arrows and a bracket indicate multiple bands of amplitude peaks of tsunami wave trains.
time window (with tapering on the 5% edges) (green line in Figure 3a).m i in Equation 1 now represents the relative source amplitude of the ith source event to that of Se02, and we impose m i ≥ 0 under our hypothesis of multiple similar events.
Using the Green's function, we deconvolve the stacked tsunami waveform (blue line in Figure 3a) to estimate the tsunami STF, following Kikuchi and Kanamori (1982).Denoting the stacked tsunami waveform as x(t), we first take a single event and determine m 1 and t 1 by minimizing the error defined as where T is the length of the stacked waveform (t = 0 at 18:30, and t = T at 23:30), obtain the residual waveform x (1) (t): (3) In the ith iteration, m i and t i are determined by minimizing the error Δ i : which yields the residual waveform x (i) (t): (5) We iteratively repeat the procedure above until the approximation error changes by less than 2% by an iteration The approximation accuracy is quantified by the normalized approximation error: Thus, we determine the relative source amplitudes of m i at the timings of t i (i = 1, …, N).
We find that, in the iterative deconvolution process, the source amplitudes determined in earlier iterations tend to be larger, because the original observed waveform is fit mainly by earlier-determined events, and later events are fit to the residual waveforms (Kikuchi & Kanamori, 1982).As seen in the iterative deconvolution results (Figure S5 in Supporting Information S1), a source event at 21:17, determined in the first iteration, has a very large source amplitude.To examine how reliable the results are, we re-determine the source amplitudes by an additional least-squares method.While fixing the source event times t i , we re-estimate the source amplitude   ′  by minimizing the following error,  () − ∑  =1  ′  ( − ) , by the non-negative least-squares method (Lawson & Hanson, 1974).Thus, we obtain the tsunami STF, represented by the relative source amplitude   ′  at the source times of t i .The tsunami waveform is modeled with Equation 1, where m i is replaced by   ′  .The additional least-squares method improves the amplitude balance of several source events close in time, determines the amplitude of an event as zero, which we remove from the event list; then, the normalized approximation error is reduced significantly from 0.174 to 0.118 (compare Figure S5 in Supporting Information S1 with Figure 3).

Results
As main results, we show the tsunami STF based on the analysis of the stacked waveform of DONET1 (Figure 3).The tsunami STF is composed of 23 source events that span from 19:53 to 22:02, labeled as Ts01-23 (Figure 3c and Table S3 in Supporting Information S1).The fit between the stacked waveform and the convolved waveform is remarkably good (Figure 3b).This supports our hypothesis that repetitive tsunami source events took place at similar locations, intermittently producing similar tsunami waveforms.
Our empirical Green's function contains about five wave cycles but excludes later waves, such as slowly propagating higher-frequency waves and reflected waves on coasts, which may cause artifacts in our estimation due to overfitting to unmodelled later waves from preceding events.Later events, such as Ts17-23, may be less reliable because the later waves from the large-amplitude preceding events (e.g., Ts13, 14 and 15) can cause nonnegligible artifacts.In Figure S6 in Supporting Information S1, we demonstrate that most parts of the stacked waveform can be well constructed by the events up to Ts16.For a similar reason, some events that are determined close in time to each other with a time difference of only <∼100 s (Ts06-07, and Ts09-10) may be separately deconvolved from a single event.
Given the imperfect Green's function, we here consider major source events took place at the timings of Ts01-05, 07-09, and 11-16, which spans from 19:53 (Ts01) to 21:26 (Ts16) for over 1.5 hr (Figure 3c).The sequence of the major events gradually increases the relative source amplitude from 1.0 to 6.5 and reduces the interval time approximately from 1,200 to 250 s (Table S3 in Supporting Information S1).
The later major events with larger amplitudes (e.g., Ts11-16) occur with inter-event times of 200-300 s (Figure 3c, and Table S3 in Supporting Information S1), which are comparable to the dominant period of the observed tsunami waveforms (Figure 2e). Figure 4a shows tsunami waveforms from Ts11-16, each of which 10.1029/2023GL106949 7 of 11 has a non-negligible amplitude of ∼10 mm, still less than half of the amplitude of the stacked largest waves.Yet, their waveform phases match with each other, and thereby the superposition of the tsunami waveforms doubles the wave amplitude and reproduces the largest waves (Figure 4b).Therefore, the late arrivals of the largest tsunami waves can be attributed to the later large events with inter-event times similar to the characteristic period of the tsunami waves.
The major events in the tsunami STF correlate well with the swarm of seismic events reported in the USGS catalog.In Figure 3d, we compare the tsunami STF with the seismic events Se01-15 (Table S1 in Supporting Information S1).Each of Se02-15, excluding Se01, nearly coincides with one of the major tsunami source events up to Ts16, and the overall trends in event size are also similar.
We also investigate origins of the T-phase signals by analyzing the vertical component of a broadband ocean bottom seismometer (OBS) at KMB06 (Figure 1b); we apply the band-pass (1-6 Hz) filter, convert the waveform into envelope by the Hilbert transform and the moving average with a 5-s window, and identify T-phase signals with a maximum amplitude larger than an empirical threshold of 2.0 × 10 3 nm/s.The origin times of the signals are estimated by shifting the maximum amplitude time backward by the T-phase travel time from a seismic source location (140.026°E,29.787°N) to KMB06 (5.48 min).As results, we detect 14 T-phase events, labeled as Tp01-14 (Figure 3d and Table S4 in Supporting Information S1).The temporal history of the T-phase events, except for Tp01, agrees with the tsunami STF (Figure 3c), as well as the seismic event swarm (Figure 3d), in terms of the origin times and the overall trend in size.
In Supporting Information S1, we show the results using the stacked tsunami waveform obtained with the 13 records of DONET2 (Figure S7 in Supporting Information S1).The estimated tsunami STF contains source events corresponding to the major events determined by the analysis of DONET1, exhibiting similar trends in size and inter-event time and the coincidence with the seismic and T-phase origin times.

Discussion and Conclusions
Table S5 in Supporting Information S1 summarizes 15 source events (labeled as EV01-15), which have similar origin times based on data of tsunami waves, seismic waves, and/or T-phases.The event EV01 at 18:58, detected as seismic and T-phase events (Se01 and Tp01), did not generate measurable tsunami waves.EV02-14, from 19:53 to 21:21, are all detected as tsunami, seismic and T-phase sources.EV15 at 21:26, listed as seismic and tsunami sources, did not excite a T-phase signal strong enough to be identified based on our detection, but its T-phase signal can be seen as a small peak in the OBS record (around Se15 in Figure 3d).Thus, we suggest that the 14 repetitive events (EV02-15) excited strong T-phases and notable tsunami waves, with seismic radiation equivalent to m b 4-5.Although the mechanism of the repetitive events remains to be solved, the coincidental excitation of T-phases and tsunamis implies a very shallow source depth in the crust or just on the seafloor.T-phase data may constrain the water depth where the events took place, considering the efficient excitation in a water depth of ∼1,000 m, a range of the so-called SOFAR channel (Okal, 2008).On the other hand, the tsunamigenesis requires a large volume of displaced seawater by such as seafloor deformation or mass movement on the seafloor.
Prompt reports have been published that indicate associations of the sequence to submarine volcanic activities.In late October 2023, Japan Coast Guard (2023) found drifted materials in the sea around Torishima and Sofugan volcanoes, parts of which were analyzed and confirmed as volcanic pumices (Earthquake Research Institute, the University of Tokyo, 2023; Geological Survey of Japan, National Institute of Advanced Industrial Science and Technology, 2023).Mizutani and Melgar (2023) analyzed parts of the tsunami waveforms, suggesting its possible origin of submarine volcanic activity.In early November, Japan Agency for Marine-Earth Science and Technology (2023) surveyed the seafloor near the source region and preliminarily reported the existence of a submarine caldera ∼20 km to the west of Sofugan volcano.Hence, we speculate the tsunami sequence is associated with volcanic unrest, most likely near the submarine caldera.
Another series of abnormal tsunamis have recurrently taken place, almost every 10 years, due to M w 5.4-5.8volcanic earthquakes at another submarine caldera, Sumisu Caldera (Figure 1c) (Kanamori et al., 1993;Satake & Kanamori, 1991), for which a submarine trapdoor faulting in the inflating caldera was proposed (Sandanbata et al., 2022); the same mechanism has been recently proposed for tsunamigenic earthquakes at two other submarine calderas: Curtis Caldera, New Zealand (Sandanbata et al., 2023), and Kita-Ioto Caldera, Japan (Sandanbata & Saito, 2024).These submarine trapdoor faulting events are tsunamigenic but exhibit atypically inefficient seismic excitation (Sandanbata et al., 2021(Sandanbata et al., , 2022)).Repetitive T-phase events possibly related to volcanic activity also followed the event in 1996 (M w 5.7) at Sumisu Caldera (Sugioka et al., 2000).Considering these past observations, the trapdoor faulting mechanism at the newly found caldera near Sofugan volcano may be a plausible candidate for the October 2023 tsunamis.However, the mechanisms proposed above cannot be distinguished from our analysis focusing on the temporal history of the tsunami generation.In future studies, it is critical to reveal the exact locations and initial shapes of the repetitive tsunamis, and a compilation of data sets of tsunamis, seismic waves, and T-phases would be the key to understand the relationship between the October 2023 sequence and submarine volcanic activity.
The repetitive character of EV02-15 infers how the series of phenomena proceeded.As seen in Figures 4c and 4d, the source event number increases exponentially with time, and the inter-event interval of the events gradually decreases.Hotovec et al. (2013) reported that the swarm of repeating earthquakes at Redout Volcano, Alaska, accelerated into an explosive eruption, with decreasing inter-event times and increasing event magnitudes.Different geological phenomena with similar characteristics were also reported for landslide events (Yamada et al., 2016), collapse events of volcanic calderas (Michon et al., 2009;T. A. Wang et al., 2023), and fault slip in tectonic environment (Igarashi, 2000).In contrast to EV02-15, EV01 did not generate measurable tsunami waves but accompanied strong T-phases and detectable seismic waves.The envelope shape of the T-phase of EV01 exhibits slower amplitude increase and longer duration than those of the others (Figure 3d), indicating its different source mechanism or properties.Although we have considered that tsunami source events Ts17-23, after EV15, are less reliable due to the imperfect Green's function, it is also possible that these later sources have generated tsunamis without notable seismic waves or T-phases, since strength of seismic waves and T-phases can significantly change due to differences in source depth and/or location (e.g., Fukao et al., 2018;Wech et al., 2018).
The 8 October 2023 tsunami sequence was preceded by another seismic swarm mainly from 2 to 6 October, including two M w > 6 earthquakes (Figure 1c).The relationship between the two swarms is unclear, but we note that some earthquakes on 5 October radiated much stronger seismic waves than the 8 October sequence, although high-frequency T-phase signals are larger for the latter (compare Figures S1 and S8 in Supporting Information S1), indicating a significant difference in their source mechanism or depth.One possible hypothesis for the link is that the preceding swarm was related to the movement of magma in the crust, leading to another phase of tsunamigenic volcanic process, such as volcanic deformation, underwater eruptions, or flank collapses.Another is that ground shaking due to the preceding swarm destabilized parts of sloped bathymetry, leading to mass failures a few days later.
The enigmatic tsunamis on 8 October 2023 sheds light on the difficulty in forecasting tsunamis without any significant earthquake in the Izu-Bonin region, where no offshore tsunami observation system is deployed.The coincidence of strong T-phases with tsunami generation may help forecast tsunamis based on the detected T-phases, as has been long suggested (Ewing et al., 1950;Matsumoto et al., 2016).However, the Izu-Bonin region hosts tens of volcano islands and submarine volcanoes, and active back-arc rift systems (Kodaira et al., 2007), implying various types of potential tsunami hazards.As discussed above, M > 5 trapdoor faulting earthquakes can cause notable tsunamis from submarine calderas (Sandanbata et al., 2022;Sandanbata & Saito, 2024).An M w 6.4 normal faulting earthquake on 24 October 2006 in a region between Sofugan volcano and Nichiyo Seamount (Figure 1c) also generated about 10-cm tsunamis (Japan Meteorological Agency, 2023).In this context, there is an urgent need to improve preparedness for tsunami hazards in the oceanic region.

Figure 1 .
Figure 1.Maps of the study area.(a) Philippine Sea off southwestern Japan region.Red inverted triangles represent locations of tide gauges with maximum tsunami heights reported by Japan Meteorological Agency (JMA).(b) Orange inverted triangles represent the Dense Oceanfloor Network (DONET) stations; each DONET node (e.g., KMA or MRA) contains 4 or 5 OBP gauges (each triangle).The tsunami is expected to arrive the earliest at KMC21 and MRE20 among DONET1 and DONET2, respectively.The ocean bottom seismometer (OBS) record at KMB06 is used for the T-phase analysis.(c) The region near Sofugan volcano.Circles represent the locations of seismic events at depths of <20 km from 1-8 October 2023, reported in the USGS earthquake catalog; open and closed circles indicate those before and after 18:00 on 8 October, respectively.Red triangles represent active volcanoes documented by Hydrographic and Oceanographic Department, Japan Coast Guard (2006).

Figure 2 .
Figure 2. Ocean bottom pressure (OBP) waveforms recorded by Dense Oceanfloor Network.(a) Waveforms after removing the tidal component by polynomial approximation, and (b) after the tidal-component removal and the band-pass filter (0.00125-0.02Hz).The amplitudes are in unit of water wave height [mm] converted from pressure, using 1.0 [Pa] = 0.102 [mmH 2 O].(c) The stacked tsunami waveform (thick black line) and the 16 waveforms of DONET1 used for the stacking (thin gray line; TableS2in Supporting Information S1).The time is aligned with the tsunami arrival at KMC21 (Figure1b).(d) The stacked waveform (average amplitude; thick black line) ± standard deviation (shaded by gray).Wavelet analysis by the continuous wavelet transformation in MATLAB(Lilly, 2017).Red arrows and a bracket indicate multiple bands of amplitude peaks of tsunami wave trains.

Figure 3 .
Figure 3.The tsunami source time function (STF) by the 24 iterative deconvolutions and the non-negative least-squares method using DONET1.(a) The empirical Green's function (green line) obtained from the stacked tsunami waveform of (blue line).(b) The convolved tsunami waveform (red line) compared with the stacked waveform.(c) Tsunami STF, composed of 23 source events (Ts01-Ts23; TableS3in Supporting Information S1).The source amplitudes are relative to that of Ts01.(d) The temporal history of the seismic (Se01-15) and T-phase events (Tp01-14), and the envelope of the vertical component of ocean bottom seismometer (OBS) at KMB06 (Figure1b), shifted backward in time assuming the T-phase speed at 1.5 km/s (gray shade).

Figure 4 .
Figure 4. Tsunami waves of the later events.(a) The waveforms due to each of Ts11-16, shown in the inset panel.(b) The convolved waveform with Ts11-16 (orange line) and the stacked waveform of DONET1 (blue line).(c, d) The relationships between (c) the event number and timing, determined by tsunamis, seismic waves, and/ or T-phases, and (d) between the event number and the inter-event time of the seismic events.