Earthquake Rupture and Tsunami Generation of the 2015 Mw 5.9 Bonin Event Revealed by In Situ Pressure Gauge Array Observations and Integrated Seismic and Tsunami Wave Simulation

On September 1, 2015, an Mw 5.9 interplate earthquake occurred near the Bonin Trench. An array of in situ ocean‐bottom absolute pressure gauges (APGs) observed its tsunami generation field consisting of static and dynamic pressure changes due to tsunami and crustal deformation and due to seismic motion, respectively, under much higher station density than ever reported. We propose an approach to synthesize the pressure change inside the focal area, which reproduces the APG waveforms well. We further successfully estimate the finite fault model of the earthquake and constrain the rupture duration only from the APG data. The relatively low stress drop and seismic wave radiation from the fault model may suggest the nature of a tsunami earthquake. The in situ APGs have a large potential to reveal the broadband source process of the earthquake, which is essential to understand the seismotectonics in the subduction zone.

considered as noises, although recent studies have utilized these signals as seafloor seismograms to estimate centroid moment tensor (CMT, An et al., 2017;Kubota et al., 2017). Furthermore, when APGs are installed inside the focal area, seismic and tsunami signals overlap, and thus it is not straightforward to utilize the pressure data as simple seismograms (Saito & Tsushima, 2016). Kubota, Saito, Tsushima et al. (2021) developed a method applicable to APGs to decompose these signals. These advancements have enabled us to utilize the broadband wave signals of the APG.
On September 1, 2015, an earthquake occurred near the Bonin Trench (hereafter, the mainshock; Figure 1a,M w 5.8,Global CMT [GCMT]; 8 km, U.S. Geological Survey [USGS]). Its reverse-faulting mechanism with a shallow dip angle indicates the interplate earthquake between the Pacific Plate and Philippine Sea Plate. In the Izu-Bonin subduction zone, the convergence rate of the Pacific Plate is ∼45 mm/year (DeMets et al., 2010), but it is thought that the two plates are weakly coupled and most slip deficit is aseismically released (Bilek & Lay, 2018;Scholz & Campos, 1995;Uyeda & Kanamori, 1979). Similar M ∼6 interplate events often occur but almost no major (M > ∼7) interplate events have been reported (e.g., Kaiho, 1991;Plescia & Hayes, 2020). Some studies point out the possibility of M ≥ 9 earthquakes in this region (Ikuta et al., 2015;McCaffrey, 1997;Plescia & Hayes, 2020). This M ∼6 mainshock will be helpful in understanding the physics of the characteristics of the seismogenesis in the Izu-Bonin subduction zone.
When the mainshock occurred, we had installed an array of APGs just above the focal area under a much higher station density than reported ever ( Figure 1a; Fukao et al., 2021), and the array recorded the tsunami generation in detail. In this study, through our integrated analysis of tsunamis and seismic waves, we estimate the mainshock rupture process and show the great advantages of the in situ APGs for seismological studies.

Theory for Tsunami Generation and Propagation
We briefly describe a theory of the tsunami generation and propagation (e.g., Saito, 2013Saito, , 2017Saito, , 2019. When earthquakes occur beneath seafloor, the seafloor vertically and horizontally displaces. Sea-surface also vertically displaces and then collapses and propagates due to gravity, resulting in tsunamis. The APG observes the hydrostatic pressure changes related to the tsunamis and the seafloor permanent displacement as where ρ 0 and g 0 are the seawater density and gravitational acceleration, and η(t) and u z (t) are the time history of the sea-surface height (tsunami) and seafloor vertical displacement, respectively (Tsushima et al., 2012). The APGs also observe the dynamic pressure changes related to the seafloor dynamic vertical motion. The dynamic pressure change is approximately given as where a z (t) is the seafloor vertical acceleration, and h 0 is the seawater depth (e.g., Filloux, 1982). This equation can be interpreted as an action-reaction force due to the vertically accelerating seafloor. Equation 2 holds when the wave frequency f is lower than the fundamental acoustic resonant frequency f 0 = c 0 /4h 0 (f < f 0 , c 0 ∼ 1.5 km/s: ocean-acoustic wave velocity), enough to consider the seawater to be incompressible (e.g., Nosov & Kolesov, 2007;Saito, 2019). Since the hydrostatic and dynamic pressure changes are generated with and without gravity, respectively, the tsunami generation field can be expressed by the sum of them (Saito, 2019;Saito & Tsushima, 2016):  (Iwasaki et al., 2015). (c, d) Comparisons of the observed (black) and simulated waveforms. Blue, green, and red traces are the synthetic waveforms considering only hydrostatic pressure change, only dynamic pressure change, and both, respectively. Low-pass filters with a cutoff of 0.033 Hz are applied. (e, f) Fourier spectra of the waveforms at B09 and B10. The black, gray, and colored lines are the spectra for the observation, background, and simulations, respectively. Colored arrows show the frequency bands where the hydrostatic and dynamic components are dominant.
The APGs recorded large pressure fluctuations of up to ∼0.05 MPa. Small pressure fluctuations, observed ∼1 min before the focal time, are due to the foreshock (the epicenter determined by Fukao et al., 2021 is shown by small white star in Figure 1a). We also show the waveforms from the broadband ocean-bottom seismometer (BBOBS) and the differential pressure gauge (DPG) at B05 in Figure S1a. The BBOBS waveform saturated, whereas the DPG successfully recorded the whole waveform without the saturation. Note that the DPG was not used in our analyses because of its lower sensitivity to the longer-period components, such as the seafloor permanent deformation.
We then applied the acausal lowpass filter with a cutoff of 0.033 Hz (red traces in Figure S1b). The cutoff was determined so that the frequency is enough lower than the fundamental acoustic resonant frequency in this region (f 0 ∼ 0.066 Hz, assuming h 0 = 5.5 km). Large impulsive signals related to the seismic wave were confirmed as well as those in the band-pass-filtered waveforms (0.01-0.033 Hz, green traces). The polarities of the first and second pulses at B10 were up and down, respectively, while the opposite feature was confirmed at B09 and other stations. Low-frequency signals attributed to the tsunamis and seafloor permanent displacements are then observed, as well as in the low-pass-filtered waveforms (≤0.01 Hz, blue traces).

Forward Simulation of APG Waveforms
We forwardly synthesize the APG waveforms to investigate the contributions of the hydrostatic and dynamic pressures. We assume a rectangular planer fault near B10 (yellow rectangle in Figure 1a). For simplicity, we assume pure reverse-faulting slip (rake = 90° and strike = 180°). Fault dip angle and depth are determined based on the plate boundary model of Iwasaki et al. (2015) and Takahashi et al. (2015) (dip = 16°, center depth = 18 km, Figure 1b). The horizontal location and dimension of the fault were determined after some trial and error (length L = 15 km, width W = 15 km). Slip amount was set as D = 13 cm so that the observed waveforms are explained (seismic moment: M 0 = 8.8 × 10 17 N m, μ = 30 GPa).
Although there are some approaches based on the coupled simulations of the seismic waves and tsunamis (e.g., Kozdon & Dunham, 2013;Lotto & Dunham, 2015;Madden et al., 2021;Maeda & Furumura, 2013;Wilson & Ma, 2021), we adopted a simpler approach. The procedure is briefly explained here; see Text S1 in Supporting Information S1 for details. The hydrostatic pressure changes (Equation 1) are synthesized by a conventional tsunami calculation method from a rectangular fault (Kajiura, 1963;Okada, 1992;Satake, 2002). The dynamic pressure changes (Equation 2) are converted from the seafloor vertical accelerations a z (t), synthesized based on the elastodynamic equation with the conventional discrete wavenumber method (e.g., Herrmann, 2013) assuming a point source at the center of the rectangular fault. In the conversion, we assume a seawater density of ρ 0 = 1.03 g/cm 3 and depths h 0 , as listed in Table S1. We finally combine them to calculate the hydrostatic and dynamic pressure changes (Equation 3) and apply the same lowpass filter as applied to the observation.
In the simulation, we assume the following source time function τ(t): where T r represents the duration of the source time function. Considering the GCMT duration of 4 s, we assumed T r = 5.0 s.
Combined waveforms of the hydrostatic and dynamic pressure changes reproduced the observation (red trace in Figure 1c). The forward simulation considering only the hydrostatic pressure changes reproduce tsunamis and offset changes due to seafloor deformation, while the dynamic pressure simulation reproduces impulsive pressure change (blue and green traces in Figure 1d). We also calculate the Fourier transform of these waveforms using the time windows during −300 and 600 s from the origin time, based on the definition of Aki and Richards (2002) (black and colored lines in Figures 1e and 1f). The background spectra calculated using the time window of −3,900 to −3,000 s from the origin time are also shown in gray lines. The hydrostatic pressure changes contribute to the signals in the frequency range of <∼0.01 Hz (compare gray and blue lines), while the components in the range of >∼0.01 Hz were reproduced by the dynamic pressure changes (gray and green lines). Both hydrostatic and dynamic components contributed to reproducing the signals between around 0.005-0.02 Hz.
We assess goodness of waveform reproductivity based on the root mean square error (RMSE): where obs i E x and syn i E x is the ith data sample of the observed and synthetic waveforms, and N is the total number of the data sample. Using all the APG stations with the time window between −300 and 600 s from the origin time, we obtain RMSE = 0.32 hPa.
When assuming larger T r , the simulated dynamic pressure changes temporally elongate and the maximum amplitudes decrease (red lines in Figures 2a and 2b), leading to larger RMSEs (green or red line in Figures 2c and 2d). On the other hand, the RMSEs considering only the hydrostatic pressure change (blue) are almost identical regardless of T r . This indicates that the dynamic pressure is sensitive to duration and that incorporating dynamic pressure provides better temporal resolution for the rupture. The duration was not well constrained for T r < ∼20 s from the RMSE values because of the cutoff of the lowpass filter (0.033 Hz), although we constrained the duration as T r ∼7 s based on another approach shown later.
In the calculation, the frequency dispersion effect of the tsunami propagation was neglected, because the simulated pressure waveforms considering the dispersion showed its effect was minor ( Figure S2). Furthermore, although the relationship p = ρ 0 g 0 η was used to convert the tsunami height to the hydrostatic pressure (Equation 1), more rigorous relationship is p = ρ 0 g 0 η/cosh(kh 0 ) (k: wavenumber, e.g., Saito, 2019). Assuming the wavelength of the tsunami λ (=2π/k) from the spatial dimension of the initial tsunami height, as ∼50 km, we obtain the factor of 1/cosh(kh 0 ) ∼ 0.8 (assuming h 0 = 5.5 km). Although we neglect this effect in the finite fault modeling shown below, this may suggest the actual seismic moment is possibly slightly larger than the estimated one.

Finite Fault Inversions
We estimated the coseismic slip distribution of the mainshock by inverting the APG data  see Text S2 in Supporting Information S1 for details). To calculate the pressure changes from each subfault (Green's function), we assume a planar fault with a dimension of 40 km length × 28 km width and divide it into 4 × 4 km subfaults. We calculate the Green's functions related to the hydrostatic pressure changes using each rectangular subfault. We also calculate the dynamic pressure change Green's functions assuming a point source at the center of each subfault. We finally combined them to calculate the Green's functions. We assume the same fault mechanism and duration as used in the forward simulation. We assume the ruptures of all subfaults begin simultaneously (i.e., infinite rupture propagation velocity), because the filter cutoff (0.033 Hz) is sufficiently low for the M ∼6 earthquake to neglect the rupture propagation across the fault. We impose a spatial smoothing constraint ( Figure S3) and a nonnegativity constraint in the inversion analysis (Lawson & Hanson, 1974).
The estimated slip distribution is plotted in Figure 3a. The location of the large slip area was consistent with the rectangular fault in the forward modeling, located at ∼15 km east and ∼10 km northwest from the USGS and GCMT centroids, respectively. The seismic moment of M 0 = 9.5 × 10 17 N m (M w 5.9) was consistent with the USGS and GCMT solutions. The simulated waveforms from the slip model (red trace in Figure 3c) reproduced the observation (RMSE = 0.30 hPa).
To investigate the rupture characteristics of this event, we calculated the stress drop at each subfault from the slip distribution (Figure 3b). The subfaults with the slip amount D is larger than 20% of the maximum slip D max (marked by black lines in Figure 3a, ∼210 km 2 ) correspond to the region where the stress drop is positive. We also calculate the slip-weighted average of the stress drop, Δσ E (Noda et al., 2013). Using the subfaults with D > 0.2D max , we obtain Δσ E = 0.5 MPa. We evaluate the estimation uncertainty using the mean and standard deviation based on the jackknife inversion test (μ ± 1σ, Figures S4 and S5, see Text S3 in Supporting Information S1 for details). We obtain the uncertainties of M 0 = 8.0-11.6 × 10 17 N m, the slip area of 208-224 km 2 , and Δσ E = 0.46-0.50 MPa.
We also conduct inversions considering only either the hydrostatic or dynamic pressure changes. When considering only the Green's function related to hydrostatic pressure changes, the slip distribution resembles the original inversion (RMSE = 0.49 hPa, Figure S6), indicating that the hydrostatic pressure contributes to the spatial constraint of the horizontal location and the horizontal extent of the rupture area. This is expected because tsunamis contain the information of the fault location and rupture extent as the tsunami source distribution Kubota et al., 2018). On the other hand, the inversion considering only the dynamic pressure changes has a spatially broader slip region, suggesting a low constraint for the rupture extent (RMSE = 0.81 hPa, Figure S7). These inversions indicate that the hydrostatic components essentially provided the effective information on the constraint of the slip distribution.
Dynamic pressure changes attributed to the foreshock (gray star in Figure 1a) show identical polarities to the mainshock (Figure 3c). We multiply the synthesized waveform due to the mainshock by 0.1 and shifting ∼1 min earlier (green line in Figure S8). This waveform fits well, indicating that the foreshock has a similar faulting mechanism to the mainshock and a seismic moment equal to 0.1 times the mainshock.

Discussion and Conclusions
The rupture time history was not resolved well even if using the dynamic pressure changes, because the rupture duration was much shorter than the filter cutoff (30 s). However, ocean-acoustic wave signals (elastic waves propagating in the seawater), ranging much higher-frequency bands, possibly give us additional information. We attempt to estimate the duration focusing on higher-frequency bands. In Figure 4, we calculate the pressure change due to the ocean-acoustic waves at the seafloor, synthesized from the isotropic stress tensor (p = −σ xx = −σ yy = −σ zz , e.g., Saito & Tsushima, 2016;Saito et al., 2019) with the wavenumber-frequency method (Herrmann, 2013). We assume a point source at the center of the rectangular fault (yellow dot in Figure 1) and use the structural model based on the Preliminary Reference Earth Model (PREM; Dziewonski & Anderson, 1981) incorporating a 6-km seawater layer ( Figure S9).
When assuming the duration T r = 7 s (Equation 4), the simulation reproduced the observed waveform at B10 well, but the simulations with longer or shorter durations did not (Figure 4a). This indicates that the rupture duration of the mainshock is ∼7 s. We note that the waveforms at the other stations were not reproduced well (Figure 4b), probably due to simple assumptions such as the point source and the structure model. To accurately reproduce the ocean-acoustic waves and thus to utilize the ocean-acoustic wave signals for the fault modeling, an appropriate structure model optimized for the target region and incorporation of the rupture propagation across the fault are necessary, even for moderate (M ∼ 6) earthquakes.
The mainshock occurred at the plate boundary near the trench axis. It is often reported that "tsunami earthquakes" (e.g., Kanamori, 1972;Tanioka & Satake, 1996) occur at the shallower portion of the plate boundary near the trench axis (Fukao, 1979). A relatively small stress drop is another feature of tsunami earthquakes (Bilek et al., 2016;Ye et al., 2016). The stress drop of the mainshock (Δσ E = 0.46-0.50 MPa) is slightly small compared with those expected in ordinary interplate earthquakes (∼10 0 MPa). Considering the fault length of ∼15 km and duration of T r ∼ 7 s, the rupture velocity across the fault can be approximated as ∼2 km/s supposing unilateral rupture propagation. We further estimate the radiation energy E R and radiation efficiency E R /M 0 (Venkataraman & Kanamori, 2004), assuming the ω −2 moment rate spectrum model (Aki, 1967) with a corner frequency of f c = 1/T r = 0.14 Hz (see Text S4 in Supporting Information S1 for details). We obtain E R = 1.77 × 10 12 N m and E R /M 0 = 1.87 × 10 −6 . When taking account of the uncertainty of these estimations, the upper limits of the possible ranges are estimated as E R ≤ 2.7 × 10 12 N m and E R / M 0 ≤ 3.3 × 10 −6 (see Text S4 in Supporting Information S1 for details). The stress drop and radiation efficiency seem slightly smaller than those for the typical interplate earthquakes but are rather consistent with those for tsunami earthquakes estimated from the global M > ∼7 event catalog (Ye et al., 2016).
Two significant aseismic slip events were also detected in and around the rupture area of the mainshock, just after and 3.5 days after the mainshock . Their magnitudes were larger than the mainshock (M w ∼ 6.5) and had time scales of ∼1 hr, which is much shorter than the characteristic duration of the typical aseismic events in this magnitude range (∼days, Ide et al., 2007). Fukao et al. (2021) interpreted these aseismic events as the transitional regime between the unstable seismic-slip and stable plate-sliding regimes (Barbot, 2019). Taking the overlapping rupture areas of these aseismic slips and the mainshock into account, the mainshock rupture might have occurred under the unstable seismic-slip regime but relatively close to the transitional regime. This might indicate that the mainshock rupture was more likely to be a tsunami earthquake (unstable but close to the transitional regime) rather than an ordinary interplate earthquake (completely unstable). The mainshock rupture characteristics may be one important feature to understand the Bonin subduction zone and may be relevant to the feature of the Bonin subduction zone such as the very small coupling rate as well as and the low normal stress along the plate boundary (Scholz & Campos, 1995).
Our analyses suggested that the pressure changes due to tsunamis have an advantage in constraining earthquake source dimensions (i.e., fault length and width) and thus the stress drop, while the use of the high-frequency pressure changes due to dynamic pressure changes and ocean-acoustic waves will give constraints on the rupture duration. Particularly, the dynamic pressure components (<∼0.03 Hz) will contribute to constraining the time history for megathrust (M ∼ 9) earthquakes which have much longer rupture duration of >10 2 s. We also point out another advantage of our approach, whereby the high-frequency tsunami signals ranging ∼10 −2 -10 −1 Hz can be used for fault modeling (Figures 1e and 1f). Conventional tsunami analyses cannot deal with tsunamis in the frequency range where low-frequency seismic wave signals are overlapped. The analysis incorporating the dynamic pressure changes makes it possible to include such high-frequency tsunami signals and estimate broadband rupture characteristics of offshore earthquakes (Webb, 1998;Webb & Nooner, 2016).
It is often difficult to constrain the slip distribution of M ∼ 6 offshore earthquakes only from seismographs or a few APGs located far from the source region (e.g., . However, our in situ APG array enabled us to reveal the detailed mainshock rupture process. Recently, new wide offshore networks, which can observe the seismicity in the wider region with uniform quality, have been established (Aoi et al., 2020;Kaneda et al., 2015;Kawaguchi et al., 2015). On the other hand, a small APG array, as used in this study, has the advantage of revealing the smaller-scaled geodynamic phenomena. APGs have a strong advantage in that their signal never saturates, whereas the BBOBS signal often saturates ( Figure S1a). Not only the wide permanent observation network, but the small temporary observation array, designed for observing specific phenomena, will be also important to understand the geodynamics in the subduction zone with high resolution. This study showed the great potential of the in situ broadband APG array data, in which the signal never saturates, to deepen our understandings of the subduction zone processes.