Confidential manuscript submitted to Journal of Geophysical Research: Solid Earth

1

Earthquake rupture and tsunami generation of the 2015 Mw 5.9 Bonin event revealed by in-situ

2pressure gauge array observations and integrated seismic and tsunami wave simulation

3
4

Tatsuya Kubota1, Tatsuhiko Saito1, Yoshio Fukao2, Hiroko Sugioka3, Aki Ito2, Takashi

5

Tonegawa2, Hajime Shiobara4, and Mikiya Yamashita5671National Research Institute for Earth Science and Disaster Resilience, Tsukuba, Japan

8

2

Japan Agency for Marine-Earth Science and Technology, Yokohama, Japan

9

3

Kobe University, Kobe, Japan

10

4

Earthqueke Research Institute, The University of Tokyo, Tokyo, Japan

11

5

Geological Survey of Japan, National Institute of Advanced Industrial Science and Technology,

12

Tsukuba, Japan

13
14

Corresponding author: T. Kubota (kubotatsu@bosai.go.jp)

15
16
17

Key Points:•

where long-period seismic waves overlap with tsunamis.

18
19

•

22
23

A numerical simulation for a fluid-elastic medium was conducted for the pressure change
synthetics above the seafloor.

20
21

We examined pressure changes inside the source region of the 2015 Bonin Mw 5.9 EQ•

Analyzing both tsunamis and dynamic ground motions enabled us to estimate the stressdrop and the rupture duration of the Mw 5.9 event.

Confidential manuscript submitted to Journal of Geophysical Research: Solid Earth24

Abstract25

On 1 September 2015, an Mw 5.9 interplate earthquake occurred near the Bonin Trench. An array26

of in-situ ocean-bottom absolute pressure gauges (APGs) observed its tsunami generation field27

consisting of static and dynamic pressure changes due to tsunami and crustal deformation and to28

seismic motion, respectively, under much higher station density than ever reported. We propose

29

an approach to synthesize the pressure change inside the focal area, which reproduces the APG

30

waveforms well. We further successfully estimated the finite fault model of the earthquake and31

constrain the rupture duration only from the APG data. The relatively low stress drop and32

seismic wave radiation from the fault model may suggest the nature of a tsunami earthquake. The

33

in-situ APGs have a large potential to reveal the broadband source process of the earthquake,

34

which is essential to understand the seismotectonics in the subduction zone.

35
36

Plain Language Summary

37

Absolute pressure gauges (APGs) installed on the seafloor have been widely used to observe

38

various phenomena such as tsunamis. When an APG is located inside a sub-seafloor earthquake

39

source region, the tsunami generation field is observed, which includes not only tsunamis but

40

also dynamic pressure changes related to seafloor dynamic motion during the earthquake fault

41

slip. This study analyzed densely distributed APG array data, which contained the tsunami

42

generation field of a magnitude-5.9 earthquake that happened near the Izu-Bonin Trench on 1

43

September 2015. This observation was performed by APG stations at higher station density than44

any data thus far reported. This study synthesized this APG array data using the theory for45

tsunami generation and propagation, and the result explained the data dramatically well. This46

indicates the validity of the tsunami generation theory for actual observation inside the47

earthquake source area. We also estimated the slip distribution across the fault of this earthquake,48

which could not previously be obtained. Our results indicate the combination of APG data49

obtained inside the focal area and the theory for tsunami generation and propagation is useful for50

extracting the detailed information on the earthquake rupture process.

51

Confidential manuscript submitted to Journal of Geophysical Research: Solid Earth52

1 Introduction

53

Recent developments of offshore seismic and geodetic observations have advanced our54

understanding of the ocean (e.g., Bürgmann & Chadwell 2014; Favali et al., 2015). The ocean-55

bottom absolute pressure gauge (APG) has also driven fundamental studies of tsunamis (e.g.,

56

Fukao et al., 2018; Levin & Nosov, 2009; Rabinovich & Eblé, 2015; Saito & Kubota, 2020) and57

seafloor geodesy (e.g., Wallace et al., 2016). Recent wide and dense offshore observation

58

networks will drive such studies more actively (Aoi et al., 2020; Kaneda et al., 2015; Kawaguchi

59

et al., 2015).

60

However, when applying conventional approaches to such new offshore data, problems

61

often happen. For example, as the APGs have been deployed above earthquake source regions, it

62turned out the seafloor permanent deformation must be considered to precisely evaluate tsunamis

63

in the APGs (Tsushima et al., 2012). APGs also observe the pressure components related to64

seismic waves (e.g., Filloux, 1982; Ito et al., 2020; Kubota et al., 2020a; Matsumoto et al., 2012;

65

Nosov & Kolesov, 2007), which have been considered as noises, although recent studies have66

utilized these signals as seafloor seismograms to estimate centroid moment tensor (CMT) (An et67

al., 2017; Kubota et al., 2017). Furthermore, when APGs are installed inside the focal area,

68

seismic and tsunami signals overlap, and thus it is not straightforward to utilize the pressure data69

as simple seismograms (Saito & Tsushima, 2016). Kubota et al. (2021) developed a method

70

applicable to APGs to decompose these signals. These advancements have enabled us to utilize

71

the broadband wave signals of the APG.

72

On 1 September 2015, an earthquake occurred near the Bonin Trench (hereafter, the73mainshock; Figure 1a, Mw 5.8, Global CMT [GCMT]; 8 km, U.S. Geological Survey [USGS]).

74

Its reverse-faulting mechanism with a shallow dip angle indicates the interplate earthquake75

between the Pacific Plate and Philippine Sea Plate. In the Izu-Bonin subduction zone, the76

convergence rate of the Pacific Plate is ~45 mm/yr (DeMets et al. 2010), but it is thought that the77

two plates are weakly coupled and most slip deficit is aseismically released (Bilek & Lay, 2018;78

Scholz & Campos, 1995; Uyeda & Kanamori, 1979). Similar M~6 interplate events often occur

79

but almost no major (M >~7) interplate events have been reported (e.g., Kaiho, 1991; Plescia &80

Hayes, 2020). Some studies point out the possibility of M ≥ 9 earthquakes in this region (Ikuta et81al., 2015; McCaffrey, 1997; Plescia & Hayes, 2020). This M~6 mainshock will be helpful inConfidential manuscript submitted to Journal of Geophysical Research: Solid Earth82

understanding the physics of the characteristics of the seismogenesis in the Izu-Bonin subduction

83

zone.

84

When the mainshock occurred, we had installed an array of APGs just above the focal

85

area under a much higher station density than reported ever (Figure 1a, Fukao et al., 2021), and86

the array recorded the tsunami generation in detail. In this study, through our integrated analysis

87

of tsunamis and seismic waves, we estimate the mainshock rupture process and show the great88

advantages of the in-situ APGs for seismological studies.

89Confidential manuscript submitted to Journal of Geophysical Research: Solid Earth90
91

Figure 1. (a) Forward simulation results. Inverted triangles denote the stations (orange: APG

92

available, open: unavailable). Yellow rectangle and CMT denote the location and mechanism of93

the fault supposed for the forward simulation, respectively. Seafloor vertical displacement94

expected from the supposed fault is shown by colors and contours (1-cm interval). The GCMT

95

and USGS solutions and the USGS epicenter are shown in gray. The small white star denotes theConfidential manuscript submitted to Journal of Geophysical Research: Solid Earth96

foreshock epicenter (Fukao et al., 2021). The black arrow denotes the plate convergence97

direction (DeMets et al., 2010). (b) Vertical cross section along line A-B in Figure 1a. The black98

line denotes the Pacific Plate surface (Iwasaki et al., 2015). (c,d) Comparisons of the observed

99(black) and simulated waveforms. Blue, green, and red traces are the synthetic waveforms

100

considering only hydrostatic pressure change, only dynamic pressure change, and both,101

respectively. Lowpass filters with a cutoff of 0.033 Hz are applied. (e,f) Fourier spectra of the102

waveforms at B09 and B10. The black, gray, and colored lines are the spectra for the103

observation, background, and simulations, respectively. Colored allows show the frequency

104

bands where the hydrostatic and dynamic components are dominant.

105
106

2 Data and method

107

2.1 Theory for tsunami generation and propagation

108

We briefly describe a theory of the tsunami generation and propagation. When109

earthquakes occur beneath seafloor, the seafloor vertically and horizontally displaces. Sea-

110

surface also vertically displaces and then collapses and propagates due to gravity, resulting in

111

tsunamis. The APG observes the hydrostatic pressure changes related to the tsunamis and the112seafloor permanent displacement as:

113
114

𝑝!"#$%&'(')* (𝑡) ≈ 𝜌+ 𝑔+ (𝜂(𝑡) − 𝑢, (𝑡),,

(1)

115
116

where ρ0 and g0 are the seawater density and gravitational acceleration, and η(t) and uz(t) are the117time history of the sea-surface height (tsunami) and seafloor vertical displacement, respectively

118(Tsushima et al., 2012). The APGs also observe the dynamic pressure changes related to the119seafloor dynamic vertical motion. The dynamic pressure change is approximately given as:

120
121

𝑝#"-(.)* (𝑡) ≈ 𝜌+ ℎ+ 𝑎, (𝑡),

(2)

122
123

where az(t) is the seafloor vertical acceleration, and h0 is the seawater depth (e.g., Filloux, 1982).

124

This equation can be interpreted as an action-reaction force due to the vertically accelerating

125

seafloor. Equation (2) holds when the wave frequency f is lower than the fundamental acoustic

126

resonant frequency f0 = c0/4h0 (f < f0, c0 ~1.5 km/s: ocean-acoustic wave velocity), enough toConfidential manuscript submitted to Journal of Geophysical Research: Solid Earth127

consider the seawater to be incompressible (e.g., Nosov & Kolesov, 2007; Saito, 2019). Since the128hydrostatic and dynamic pressure changes are generated with and without gravity, respectively,

129

the tsunami generation field can be expressed by the sum of them (Saito & Tsushima, 2016;130

Saito, 2019):

131

𝑝(𝑡) = 𝑝!"#$%&'(')* (𝑡) + 𝑝#"-(.)* (𝑡).

132

(3)

133
134

2.2 APG array observation

135

We installed an array of stations B01 to B10, inward of the northern Bonin Trench

136(Figure 1a, Table S1, Fukao et al., 2021). We use the APGs at B01–B03 and B06–B10 for which

137

data were available. Ocean tides were removed by fitting a set of harmonic oscillation functions

138(Figure S1a, see Fukao et al., 2021). The APGs recorded large pressure fluctuations of up to

139~0.05 MPa. Small pressure fluctuations, observed ~1 min before the focal time, are due to the140

foreshock (the epicenter determined by Fukao et al. (2021) is shown by small white star in

141

Figure 1a). We also show the waveforms from the broadband ocean-bottom seismometer

142(BBOBS) and the differential pressure gauge (DPG) at B05 in Figure S1a. The BBOBS

143

waveform saturated, whereas the DPG successfully recorded the whole waveform without the144

saturation. Note that the DPG was not used in our analyses because of its lower sensitivity to the145

longer-period components, such as the seafloor permanent deformation.

146

We then applied the acausal lowpass filter with a cutoff of 0.033 Hz (red traces in Figure147

S1b). The cutoff was determined so that the frequency is enough lower than the fundamental

148

acoustic resonant frequency in this region (f0 ~0.066 Hz, assuming h0 = 5.5 km). Large impulsive

149

signals related to the seismic wave were confirmed as well as those in the bandpass-filtered150

waveforms (0.01 to 0.033 Hz, green traces). The polarities of the first and second pulses at B10151

were up and down, respectively, while the opposite feature was confirmed at B09 and other152

stations. Low-frequency signals attributed to the tsunamis and seafloor permanent displacements

153

are then observed, as well as in the lowpass-filtered waveforms (≤ 0.01 Hz, blue traces).

154
155

3 Forward simulation of APG waveforms

156

We forwardly synthesize the APG waveforms to investigate the contributions of the157

hydrostatic and dynamic pressures. We assume a rectangular planer fault near B10 (yellow

Confidential manuscript submitted to Journal of Geophysical Research: Solid Earth

158

rectangle in Figure 1a). For simplicity, we assume pure reverse-faulting slip (rake = 90° and159

strike = 180°). Fault dip angle and depth are determined based on the plate boundary model of160

Iwasaki et al. (2015) and Takahashi et al. (2015) (dip = 16°, center depth = 18 km, Figure 1b).

161

The horizontal location and dimension of the fault were determined after some trial and error

162(length L = 15 km, width W = 15 km). Slip amount was set as D = 13 cm so that the observed

163

waveforms are explained (seismic moment: M0 = 8.8 × 1017 Nm, μ = 30 GPa).

164

Although there are some approaches based on the coupled simulations of the seismic165

waves and tsunamis (e.g., Kozdon & Dunham, 2013; Lotto & Dunham, 2015; Madden et al.,1662021; Maeda & Furumura, 2013; Wilson & Ma, 2021), we adopted a simpler approach. The

167

procedure is briefly explained here; see Text S1 for details. The hydrostatic pressure changes

168(Eq. (1)) are synthesized by a conventional tsunami calculation method from a rectangular fault

169(Kajiura, 1963; Okada, 1992; Satake, 2002). The dynamic pressure changes (Eq. (2)) are

170

converted from the seafloor vertical accelerations az(t), synthesized based on the elastodynamic

171

equation with the conventional discrete wavenumber method (e.g., Herrmann, 2013) assuming a172

point source at the center of the rectangular fault. In the conversion, we assume a seawater173

density of ρ0 = 1.03 g/cm3 and depths h0, as listed in Table S1. We finally combine them to174

calculate the hydrostatic and dynamic pressure changes (Eq. (3)) and apply the same lowpass

175

filter as applied to the observation.

176

In the simulation, we assume the following source time function τ(t):

177
/

12'

!

!

𝜏(𝑡) = 0 21 − cos 7 0 89,

178

(4)

179
180

where Tr represents the duration of the source time function. Considering the GCMT duration of181

4 s, we assumed Tr = 5.0 s.

182

Combined waveforms of the hydrostatic and dynamic pressure changes reproduced the183

observation (red trace in Figure 1c). The forward simulation considering only the hydrostatic

184

pressure changes reproduce tsunamis and offset changes due to seafloor deformation, while the185

dynamic pressure simulation reproduces impulsive pressure change (blue and green traces in186

Figure 1d). We also calculate the Fourier transform of these waveforms using the time windows187

during −300 and 600 s from the origin time, based on the definition of Aki and Richards (2002)

188(black and colored lines in Figures 1e and 1f). The background spectra calculated using the timeConfidential manuscript submitted to Journal of Geophysical Research: Solid Earth189

window of −3,900 to −3,000 s from the origin time are also shown in gray lines. The hydrostatic

190

pressure changes contribute to the signals in the frequency range of < ~0.01 Hz (compare gray191

and blue lines), while the components in the range of > ~0.01 Hz were reproduced by the192

dynamic pressure changes (gray and green lines). Both hydrostatic and dynamic components

193

contributed to reproducing the signals between around 0.005–0.02 Hz.

194
195

We assess goodness of waveform reproductivity based on the root mean square error(RMSE):

196
/

678 1

456
𝑅𝑀𝑆𝐸 = >3 ∑3− 𝑥))9/(𝑥)

197

, ,

(6)

198
678

199

where 𝑥)456 and 𝑥)

200

the total number of the data sample. Using all the APG stations with the time window between

201−300 to 600 s from the origin time, we obtain RMSE = 0.32 hPa.

is the i-th data sample of the observed and synthetic waveforms, and N is

202

When assuming larger Tr, the simulated dynamic pressure changes temporally elongate

203

and the maximum amplitudes decrease (red lines in Figures 2a and 2b), leading to larger RMSEs

204(green or red line in Figures 2c and 2d). On the other hand, the RMSEs considering only the205

hydrostatic pressure change (blue) are almost identical regardless of Tr. This indicates that the

206

dynamic pressure is sensitive to duration and that incorporating dynamic pressure provides better207

temporal resolution for the rupture. The duration was not well constrained for Tr < ~20 s from the208

RMSE values because of the cutoff of the lowpass filter (0.033 Hz), although we constrained the209

duration as Tr ~7 s based on another approach shown later.

210

In the calculation, the frequency dispersion effect of the tsunami propagation was

211

neglected, because the simulated pressure waveforms considering the dispersion showed its

212

effect was minor (Figure S2). Furthermore, although the relationship p = ρ0g0η was used to213

convert the tsunami height to the hydrostatic pressure (Eq. (1)), more rigorous relationship is p =214

ρ0g0η/cosh(kh0) (k: wavenumber, e.g., Saito, 2019). Assuming the wavelength of the tsunami λ (=2152π/k) from the spatial dimension of the initial tsunami height, as ~50 km, we obtain the factor of

216

1/cosh(kh0) ~0.8 (assuming h0 = 5.5 km). Although we neglect this effect in the finite fault

217

modeling shown below, this may suggest the actual seismic moment is possibly slightly larger

218

than the estimated one.

Confidential manuscript submitted to Journal of Geophysical Research: Solid Earth

219

220
221

Figure 2. Relationship between the duration and the synthesized waveforms. (a,b) Comparisons222of observed (black) and synthesized (colored) waveforms with different durations for (a) B09

223

and (b) B10. (c,d) Relationship between Tr and RMSE between the observed and simulated

224

waveforms at (c) B09 and (d) B10. Blue, green, and red traces are the synthesized waveforms

225

considering only hydrostatic, only dynamic, and both pressure changes, respectively.

226
227
228
229

4 Finite fault inversions
We estimated the coseismic slip distribution of the mainshock by inverting the APG data
(Fukao et al., 2021; see Text S2 for details). To calculate the pressure changes from eachConfidential manuscript submitted to Journal of Geophysical Research: Solid Earth230

subfault (Green’s function), we assume a planar fault with a dimension of 40 km length × 28 km

231

width and divide it into 4 km × 4 km sub-faults. We calculate the Green’s functions related to the232

hydrostatic pressure changes using each rectangular subfault. We also calculate the dynamic

233

pressure change Green’s functions assuming a point source at the center of each subfault. We

234

finally combined them to calculate the Green’s functions. We assume the same fault mechanism

235

and duration as used in the forward simulation. We assume the ruptures of all subfaults begin

236

simultaneously (i.e., infinite rupture propagation velocity), because the filter cutoff (0.033 Hz) is237

sufficiently low for the M~6 earthquake to neglect the rupture propagation across the fault. We

238

impose a spatial smoothing constraint (Figure S3) and a non-negativity constraint in the239

inversion analysis (Lawson & Hanson, 1974).

240

The estimated slip distribution is plotted in Figure 3a. The location of the large slip area241

was consistent with the rectangular fault in the forward modeling, located at ~15 km east and242~10 km northwest from the USGS and GCMT centroids, respectively. The seismic moment of243

M0 = 9.5 × 1017 Nm (Mw 5.9) was consistent with the USGS and GCMT solutions. The simulated

244

waveforms from the slip model (red trace in Figure 3c) reproduced the observation (RMSE =245

0.30 hPa).

246

To investigate the rupture characteristics of this event, we calculated the stress drop at

247

each subfault from the slip distribution (Figure 3b). The subfaults with the slip amount D is248

larger than 20% of the maximum slip Dmax (marked by black lines in Figure 3a, ~210 km2)

249

correspond to the region where the stress drop is positive. We also calculate the slip-weighted

250

average of the stress drop, ΔσE (Noda et al., 2013). Using the sub-faults with D > 0.2Dmax, we

251

obtain ΔσE = 0.5 MPa. We evaluate the estimation uncertainty using the mean and standard

252

deviation based on the jackknife inversion test (μ ± 1σ, Figures S4 and S5, see Text S3 for253

details). we obtain the uncertainties of M0 = 8.0–11.6 × 1017 Nm, the slip area of 208–224 km2,

254

and ΔσE = 0.46–0.50 MPa.

255

Confidential manuscript submitted to Journal of Geophysical Research: Solid Earth

256
257

Figure 3. Result of the inversion. (a) Slip distribution (colored tiles) and seafloor vertical258

displacement (black contours). The area of D > 0.2Dmax is marked by thick black lines. The

259

yellow circle and CMT solution denote the location of the centroid calculated from the slip

260

distribution and the assumed fault mechanism, respectively. The USGS and GCMT centroids are

261

shown by gray circles. (b) Stress drop distribution. (c) Comparisons of the inverted (black) and262

synthesized (red) waveforms. (d) Decomposition of the forward-calculated waveforms (blue:

263

hydrostatic, green: dynamic pressure).

264
265

We also conduct inversions considering only either the hydrostatic or dynamic pressure

266changes. When considering only the Green’s function related to hydrostatic pressure changes, the267

slip distribution resembles the original inversion (RMSE = 0.49 hPa, Figure S6), indicating that

Confidential manuscript submitted to Journal of Geophysical Research: Solid Earth268

the hydrostatic pressure contributes to the spatial constraint of the horizontal location and the269

horizontal extent of the rupture area. This is expected because tsunamis contain the information

270of the fault location and rupture extent as the tsunami source distribution (An et al., 2018;

271

Kubota et al., 2018). On the other hand, the inversion considering only the dynamic pressure272changes has a spatially broader slip region, suggesting a low constraint for the rupture extent

273(RMSE = 0.81 hPa, Figure S7). These inversions indicate that the hydrostatic components

274

essentially provided the effective information on the constraint of the slip distribution.

275

Dynamic pressure changes attributed to the foreshock (gray star in Figure 1a) show

276

identical polarities to the mainshock (Figure 3c). We multiply the synthesized waveform due to277

the mainshock by 0.1 and shifting ~1 min earlier (green line in Figure S8). This waveform fits278

well, indicating that the foreshock has a similar faulting mechanism to the mainshock and a279

seismic moment equal to 0.1 times the mainshock.

280
281
282

5 Discussion and conclusions
The rupture time history was not resolved well even if using the dynamic pressure283

changes, because the rupture duration was much shorter than the filter cutoff (30 s). However,284

ocean-acoustic wave signals (elastic waves propagating in the seawater), ranging much higher-

285

frequency bands, possibly give us additional information. We attempt to estimate the duration

286

focusing on higher frequency bands. In Figure 4, we calculate the pressure change due to the287

ocean-acoustic waves at the seafloor, synthesized from the isotropic stress tensor (p = −σxx = −σyy288

= −σzz, e.g., Saito & Tsushima, 2016; Saito et al., 2019) with the wavenumber-frequency method

289(Herrmann, 2013). We assume a point source at the center of the rectangular fault (yellow dot in290

Figure 1) and use the structural model based on the Preliminary Reference Earth Model (PREM,

291

Dziewonski & Anderson, 1981) incorporating a 6-km seawater layer (Figure S9).

292

When assuming the duration Tr = 7 s (Eq. (4)), the simulation reproduced the observed

293

waveform at B10 well, but the simulations with longer or shorter durations did not (Figure 4a).

294This indicates that the rupture duration of the mainshock is ~7 s. We note that the waveforms at295

the other stations were not reproduced well (Figure 4b), probably due to simple assumptions

296

such as the point source and the structure model. To accurately reproduce the ocean-acoustic

297

waves and thus to utilize the ocean-acoustic wave signals for the fault modeling, an appropriate

Confidential manuscript submitted to Journal of Geophysical Research: Solid Earth298

structure model optimized for the target region and incorporation of the rupture propagation299

across the fault are necessary, even for moderate (M~6) earthquakes.

300

301
302

Figure 4. Comparison of the observed and simulated waveforms of the ocean-acoustic waves.

303(b) Observed (black) and simulated waveforms (red) at B10 with different durations Tr. The

304

cross-correlation coefficients between the observation and simulation (CC) are also shown. (b)305

Waveforms at representative stations. Lowpass filters with a cutoff of 0.5 Hz are applied. The

306

origin time of the simulations is manually shifted.

307
308

The mainshock occurred at the plate boundary near the trench axis. It is often reported

309

that “tsunami earthquakes” (e.g., Kanamori, 1972; Tanioka & Satake, 1996) occur at the310

shallower portion of the plate boundary near the trench axis (Fukao, 1979). A relatively small

311

stress drop is another feature of tsunami earthquakes (Bilek et al., 2016; Ye et al., 2016). The

312

stress drop of the mainshock (ΔσE= 0.46–0.50 MPa) is slightly small compared with those

313

expected in ordinary interplate earthquakes (~100 MPa). Considering the fault length of ~15 km

314

and duration of Tr ~ 7 s, the rupture velocity across the fault can be approximated as ~2 km/s

315

supposing unilateral rupture propagation. We further estimate the radiation energy ER and

316

radiation efficiency ER/M0 (Venkataraman & Kanamori, 2004), assuming the ω−2 moment rate

317

spectrum model (Aki, 1967) with a corner frequency of fc = 1/Tr = 0.14 Hz (see Text S3 forConfidential manuscript submitted to Journal of Geophysical Research: Solid Earth318

details). We obtain ER = 1.77 × 1012 Nm and ER/M0 = 1.87 × 10−6. When taking account of the319

uncertainty of these estimations, the upper limits of the possible ranges are estimated as ER ≤ 2.7

320

×1012 Nm and ER/M0 ≤ 3.3 ×10−6 (see Text S4 for details). The stress drop and radiation321

efficiency seem slightly smaller than those for the typical interplate earthquakes, but are rather

322

consistent with those for tsunami earthquakes estimated from the global M > ~7 event catalog

323(Ye et al., 2016).

324

Two significant aseismic slip events were also detected in and around the rupture area of325

the mainshock, just after and 3.5 days after the mainshock (Fukao et al. 2021). Their magnitudes

326

were larger than the mainshock (Mw ~6.5) and had time scales of ~1 h, which is much shorter

327

than the characteristic duration of the typical aseismic events in this magnitude range (~days, Ide328

et al., 2007). Fukao et al. (2021) interpreted these aseismic events as the transitional regime

329

between the unstable seismic-slip and stable plate-sliding regimes (Barbot, 2019). Taking the330

overlapping rupture areas of these aseismic slips and the mainshock into account, the mainshock331

rupture might have occurred under the unstable seismic-slip regime but relatively close to the332transitional regime. This might indicate that the mainshock rupture was more likely to be a333

tsunami earthquake (unstable but close to the transitional regime) rather than an ordinary

334

interplate earthquake (completely unstable). The mainshock rupture characteristics may be one

335

important feature to understand the Bonin subduction zone and may be relevant to the feature of336

the Bonin subduction zone such as the very small coupling rate as well as and the low normal

337

stress along the plate boundary (Scholz & Campos, 1995).

338

Our analyses suggested that the pressure changes due to tsunamis have an advantage in339

constraining earthquake source dimensions (i.e., fault length and width) and thus the stress drop,

340

while the use of the high-frequency pressure changes due to dynamic pressure changes and341

ocean-acoustic waves will give constraints on the rupture duration. Particularly, the dynamic

342

pressure components (< ~0.03 Hz) will contribute to constraining the time history for megathrust343(M~9) earthquakes which have much longer rupture duration of > 102 s. We also point out

344

another advantage of our approach, whereby the high-frequency tsunami signals ranging

345

approximately 10-2 to 10-1 Hz can be used for fault modeling (Figures 1e and 1f). Conventional346

tsunami analyses cannot deal with tsunamis in the frequency range where low-frequency seismic347

wave signals are overlapped. The analysis incorporating the dynamic pressure changes makes itConfidential manuscript submitted to Journal of Geophysical Research: Solid Earth348

possible to include such high-frequency tsunami signals and estimate broadband rupture349

characteristics of offshore earthquakes (Webb, 1998; Webb & Nooner, 2016).

350

It is often difficult to constrain the slip distribution of M~6 offshore earthquakes only

351

from seismographs or a few APGs located far from the source region (e.g., Kubota et al., 2020b).

352

However, our in-situ APG array enabled us to reveal the detailed mainshock rupture process.

353

Recently, new wide offshore networks, which can observe the seismicity in the wider region with354

uniform quality, have been established (Aoi et al., 2020; Kaneda et al., 2015; Kawaguchi et al.,3552015). On the other hand, a small APG array, as used in this study, has the advantage of356

revealing the smaller-scaled geodynamic phenomena. APGs have a strong advantage in that their357

signal never saturates, whereas the BBOBS signal often saturates (Figure S1a). Not only the358wide permanent observation network, but the small temporary observation array, designed for359

observing specific phenomena, will be also important to understand the geodynamics in the360

subduction zone with high resolution. This study showed the great potential of the in-situ

361

broadband APG array data, in which the signal never saturates, to deepen our understandings of362

the subduction zone processes.

363
364
365

Acknowledgments, Samples, and Data
We used 500 m gridded bathymetry data published by the Japan Oceanographic Data

366

Center, Japan Coast Guard (https://www.jodc.go.jp/jodcweb/JDOSS/index.html, the information367

of the data is available from https://www.jodc.go.jp/jodcweb/JDOSS/infoJEGG.html). The plate368

boundary model by Iwasaki et al. (2015) was downloaded from http://evrrss.eri.u-

369

tokyo.ac.jp/database/PLATEmodel/. We used the Computer Programs in Seismology (Herrmann,

370

2013) version 3.30 for the seismic wave propagation simulation

371(http://www.eas.slu.edu/eqc/eqccps.html). We used Seismic Analysis Code software version

372

101.6a for the data processing (Goldstein et al., 2003). We used Generic Mapping Tools version373

6.2.0 (Wessel et al. 2019) to generate the figures. The station locations of the APGs are listed in

374

Table S1 and are presented in Fukao et al. (2021). The APG data is available on the Zenodo

375

repository (https://doi.org/10.5281/zenodo.5139450), which is the supplementary dataset of376

Fukao et al. (2021).

377
378

This work was supported by the Sasakawa Scientific Research Grant #2019-2037 fromthe Japan Science Society, by JSPS KAKENHI Grant Numbers JP25247074, JP17K05646,

Confidential manuscript submitted to Journal of Geophysical Research: Solid Earth379

JP19K14818, and JP20K04142, and by ERI JURP, 2014-G-02 and 2015-G-07. We would like to380

thank T. Oi for his help with the instrumentation. We would also like to thank the shipboard381

scientists and technical staff of cruise YK15-08 of R/V Yokosuka and cruise MR16-E02 of R/V

382

Mirai. The English in this manuscript was edited by Forte Science Communications (www.forte-383

science.co.jp).

384Confidential manuscript submitted to Journal of Geophysical Research: Solid Earth

385

References

386

Aki, K. (1967). Scaling law of seismic spectrum. Journal of Geophysical Research, 72, 1217–

387
388
389
390

1231. https://doi.org/10.1029/JZ072i004p01217
Aki, K., & Richards, P. G. (2002). Quantitative seismology (2nd ed.). Mill Valley, CA:
University Science Books.
An, C., Cai, C., Zheng, Y., Meng, L., & Liu, P. (2017). Theoretical solution and applications of391

ocean bottom pressure induced by seismic seafloor motion. Geophysical Research

392

Letters, 44, 10272–10281. https://doi.org/10.1002/2017GL075137

393

An, C., Liu, H., Ren, Z., & Yuan, Y. (2018). Prediction of tsunami waves by uniform slip

394models. Journal of Geophysical Research: Oceans, 123, 8366–8382.

395

https://doi.org/10.1029/2018JC014363

396

Aoi, S., Asano, Y., Kunugi, T., Kimura, T., Uehira, K., Takahashi, N., & Ueda, H. (2020).

397

MOWLAS : NIED observation network for earthquake, tsunami and volcano. Earth,

398

Planets and Space, 72, 126. https://doi.org/10.1186/s40623-020-01250-x

399

Barbot, S. (2019). Slow-slip, slow earthquakes, period-two cycles, full and partial ruptures, and400

deterministic chaos in a single asperity fault. Tectonophysics, 768, 228171.

401

https://doi.org/10.1016/j.tecto.2019.228171

402
403
404

Bilek, S. L., & Lay, T. (2018). Subduction zone megathrust earthquakes. Geosphere, 14, 1468–
1500. https://doi.org/10.1130/GES01608.1
Bilek, S. L., Rotman, H. M. M., & Phillips, W. S. (2016). Low stress drop earthquakes in the405

rupture zone of the 1992 Nicaragua tsunami earthquake. Geophysical Research Letters,

406

43, 10180–10188. https://doi.org/10.1002/2016GL070409

407

Bürgmann, R., & Chadwell, D. (2014). Seafloor geodesy. Annual Review of Earth and Planetary

408

Sciences, 42, 509–534. https://doi.org/10.1146/annurev-earth-060313-054953

409

DeMets, C., Gordon, R. G., & Argus, D. F. (2010). Geologically current plate motions.

410

Geophysical Journal International, 181, 1–80. https://doi.org/10.1111/j.1365-

411

246X.2009.04491.x

412

Dziewonski, A. M., & Anderson, D. L. (1981). Preliminary reference Earth model. Physics of the

413

Earth and Planetary Interiors, 25, 297–356. https://doi.org/10.1016/0031-

414

9201(81)90046-7

Confidential manuscript submitted to Journal of Geophysical Research: Solid Earth415

Favali, P., Beranzoli, L., De Santis, A. (Eds.). (2015). Seafloor Observatories: a new vision of416

the Earth from the Abyss. Berlin, Heidelberg: Springer. https://doi.org/ 10.1007/978-3-

417

642-11374-1

418
419
420

Filloux, J. H. (1982). Tsunami recorded on the open ocean floor. Geophysical Research Letters,
9(1), 25–28. https://doi.org/10.1029/GL009i001p00025
Fukao, Y. (1979). Tsunami earthquakes and subduction processes near deep-sea trenches.

421

Journal of Geophysical Research: Solid Earth, 84(B5), 2303–2314.

422

https://doi.org/10.1029/JB084iB05p02303

423

Fukao, Y., Sandanbata, O., Sugioka, H., Ito, A., Shiobara, H., Watada, S., & Satake, K. (2018).

424

Mechanism of the 2015 volcanic tsunami earthquake near Torishima, Japan. Science

425

Advances, 4(4), eaao0219. https://doi.org/10.1126/sciadv.aao0219

426

Fukao, Y., Kubota, T., Sugioka, Y., Ito, A., Tonegawa, T., Shiobara, H., et al. (2021) Strain

427

release by “rapid” aseismic slip at the Izu-Bonin trench. Journal of Geophysical428

Research, 126, e2021JB022132. https://doi.org/10.1029/2021JB022132

429

Goldstein, P., Dodge, D., Firpo, M., & Minner L. (2003). SAC2000: Signal processing and430

analysis tools for seismologists and engineers. In: W. H. K. Lee, H. Kanamori, P. C.

431

Jennings, & C. Kisslinger (Eds.), International Handbook of Earthquake and432

Engineering Seismology (Vol. 81(B), pp. 1613–1614). London: Academic Press.

433

https://doi.org/10.1016/S0074-6142(03)80284-X

434

Herrmann, R. B. (2013). Computer Programs in Seismology: An Evolving Tool for Instruction

435

and Research. Seismological Research Letters, 84, 1081–1088.

436

https://doi.org/10.1785/0220110096

437
438
439

Ide, S., Beroza, G. C., Shelly, D. R., & Uchide, T. (2007). A scaling law for slow earthquakes.
Nature, 447, 76–79. https://doi.org/10.1038/nature05780
Ikuta, R., Mitsui, Y., Kurokawa, Y., & Ando, M. (2015). Evaluation of strain accumulation in440

global subduction zones from seismicity data. Earth, Planets and Space, 67, 192.

441

https://doi.org/10.1186/s40623-015-0361-5

442

Ito, Y., Webb, S. C., Kaneko, Y., Wallace, L. M., & Hino, R. (2020). Sea surface gravity waves

443

excited by dynamic ground motions from large regional earthquakes. Seismological444

Research Letters. https://doi.org/10.1785/0220190267

Confidential manuscript submitted to Journal of Geophysical Research: Solid Earth445

Iwasaki, T., Sato, H., Shinohara, M., Ishiyama, T., & Hashima, A. (2015). Fundamental structure

446

model of island arcs and subducted plates in and around Japan. Abstract T31B‐2878

447

presented at 2015 AGU Fall Meeting, American Geophysical Union, San Francisco,448

CA, 14–18, Dec, 2015.

449

Kaiho, Y. (1991). Seismicity of the Izu-Ogasawara arc region. Journal of Geography (Chigaku

450

Zasshi), 100(4), 503–513 (in Japanese with English abstract).

451

https://doi.org/10.5026/jgeography.100.503

452
453
454
455

Kajiura, K. (1963). The leading wave of a tsunami. Bulletin of the Earthquake ResearchInstitute, 41, 535–571.
Kanamori, H. (1972). Mechanism of tsunami earthquakes. Physics of the Earth and Planetary
Interiors, 6, 346–359. https://doi.org/10.1016/0031-9201(72)90058-1

456

Kaneda, Y., Kawaguchi, K., Araki, E., Matsumoto, H., Nakamura, T., Kamiya, S., et al. (2015).

457

Development and application of an advanced ocean floor network system for458

megathrust earthquakes and tsunamis. In P. Favali, L. Beranzoli, & A. De Santis (Eds.),

459

Seafloor Observatories: A new vision of the Earth from the Abyss (pp. 643–662). Berlin,

460

Heidelberg: Springer. https://doi.org/10.1007/978-3-642-11374-1_25

461

Kawaguchi, K., Kaneko, S., Nishida, T., & Komine, T. (2015). Construction of the DONET real-

462

time seafloor observatory for earthquakes and tsunami monitoring. In P. Favali, L.

463

Beranzoli, & A. De Santis (Eds.), Seafloor Observatories: A New Vision of the Earth464

from the Abyss (pp. 211–228). https://doi.org/10.1007/978-3-642-11374-1_10

465

Kozdon, J. E., & Dunham, E. M. (2013). Rupture to the trench: dynamic rupture simulations of466

the 11 march 2011 Tohoku earthquake. Bulletin of the Seismological Society of

467

America, 103(2B), 1275–1289. https://doi.org/10.1785/0120120136

468

Kubota, T., Saito, T., Suzuki, W., & Hino, R. (2017). Estimation of seismic centroid moment469

tensor using ocean bottom pressure gauges as seismometers. Geophysical Research

470

Letters, 44, 10907–10915. https://doi.org/10.1002/2017GL075386

471

Kubota, T., Saito, T., Ito, Y., Kaneko, Y., Wallace, L. M., Suzuki, S., et al. (2018). Using

472

tsunami waves reflected at the coast to improve offshore earthquake source parameters:

473

application to the 2016 Mw 7.1 Te Araroa earthquake, New Zealand. Journal of474

Geophysical Research: Solid Earth, 123, 8767–8779.

475

https://doi.org/10.1029/2018JB015832

Confidential manuscript submitted to Journal of Geophysical Research: Solid Earth

476

Kubota, T., Saito, T., Chikasada, N.Y., & Suzuki, W. (2020a). Ultra-broadband seismic and477

tsunami wave observation of high-sampling ocean-bottom pressure gauge covering

478

periods from seconds to hours. Earth and Space Science, 7, e2020ea001197.

479

https://doi.org/10.1029/2020EA001197

480

Kubota, T., Saito, T., & Suzuki, W. (2020b). Millimeter-scale tsunami detected by a wide and481

dense observation array in the deep ocean: fault modeling of an Mw 6.0 interplate482

earthquake off Sanriku, NE Japan. Geophysical Research Letters, 47, e2019GL085842.

483

https://doi.org/10.1029/2019GL085842

484

Kubota, T., Saito, T., Tsushima, H., Hino, R., Ohta, Y., Suzuki S, & Inazu D. (2021). Extracting

485

near-field seismograms from ocean-bottom pressure gauge inside the focal area:

486

application to the 2011 Mw 9.1 Tohoku-Oki earthquake. Geophysical Research Letters,

487

48, e2020GL091664. https://doi.org/10.1029/2020GL091664

488
489
490
491
492

Lawson, C. L., & Hanson, B. J. (1974). Solving least squares problems. Englewood Cliffs, NJ:
Prentice-Hall.
Levin, B. W. & Nosov, M. A. (2009). Physics of tsunamis. Netherlands: Springer Netherlands.
https://doi.org/10.1007/978-1-4020-8856-8
Lotto, G. C., & Dunham, E. M. (2015). High-order finite difference modeling of tsunami

493

generation in a compressible ocean from offshore earthquakes. Computational494

Geosciences, 19, 327–340. https://doi.org/10.1007/s10596-015-9472-0

495Madden, E. H., Bader, M., Behrens, J., Van Dinther, Y., Gabriel, A. A., Rannabauer, L., Ulrich,

496

T., Uphoff, C., Vater, S., & Van Zelst, I. (2021). Linked 3-D modelling of megathrust

497

earthquake-tsunami events: From subduction to tsunami run up. Geophysical Journal

498

International, 224, 487–516. https://doi.org/10.1093/gji/ggaa484

499

Maeda, T., & Furumura, T. (2013). FDM Simulation of Seismic Waves, Ocean Acoustic Waves,

500

and Tsunamis Based on Tsunami-Coupled Equations of Motion. Pure and Applied

501

Geophysics, 170(1–2), 109–127. https://doi.org/10.1007/s00024-011-0430-z

502

Matsumoto, H., Inoue, S., & Ohmachi, T. (2012). Dynamic response of bottom water pressure503

due to the 2011 Tohoku earthquake. Journal of Disaster Research, 7, 468–475.

504

https://doi.org/10.20965/jdr.2012.p0468

Confidential manuscript submitted to Journal of Geophysical Research: Solid Earth

505

McCaffrey, R. (1997). Influences of recurrence times and fault zone temperatures on the age-rate506

dependence of subduction zone seismicity. Journal of Geophysical Research, 102(B10),50722839–22854. https://doi.org/10.1029/97jb01827

508

Noda, H., Lapusta, N., & Kanamori, H. (2013). Comparison of average stress drop measures for509

ruptures with heterogeneous stress change and implications for earthquake physics.

510

Geophysical Journal International, 193, 1691–1712. https://doi.org/10.1093/gji/ggt074

511

Nosov, M. A., & Kolesov, S. V. (2007). Elastic oscillations of water column in the 2003

512

Tokachi-oki tsunami source: in-situ measurements and 3-D numerical modelling.

513

Natural Hazards and Earth System Science, 7, 243–249. https://doi.org/10.5194/nhess-514

7-243-2007

515
516
517

Okada, Y. (1992). Internal deformation due to shear and tensile faults in a half-space. Bulletin of
the Seismological Society of America, 82, 1018–1040.
Plescia, S. M., & Hayes, G. P. (2020). Geometric controls on megathrust earthquakes.

518

Geophysical Journal International, 222, 1270–1282.

519

https://doi.org/10.1093/gji/ggaa254

520

Rabinovich, A. B., & Eblé, M. C. (2015). Deep-Ocean Measurements of Tsunami Waves. Pure

521

and Applied Geophysics, 172, 3281–3312. https://doi.org/10.1007/s00024-015-1058-1

522

Saito, T. (2013). Dynamic tsunami generation due to sea-bottom deformation: Analytical523

representation based on linear potential theory. Earth, Planets and Space, 65, 1411–

524

1423. https://doi.org/10.5047/eps.2013.07.004

525
526
527
528

Saito, T. (2017). Tsunami generation: validity and limitations of conventional theories.
Geophysical Journal International, 210, 1888–1900. https://doi.org/10.1093/gji/ggx275
Saito, T. (2019). Tsunami generation and propagation. Tokyo: Springer Tokyo.
https://doi.org/10.1007/978-4-431-56850-6

529

Saito, T., & Tsushima, H. (2016). Synthesizing ocean bottom pressure records including seismic

530

wave and tsunami contributions: Toward realistic tests of monitoring systems. Journal531

of Geophysical Research: Solid Earth, 121, 8175–8195.

532

https://doi.org/10.1002/2016JB013195

533

Saito, T., & Kubota, T. (2020). Tsunami modeling for the deep sea and inside focal areas.

534

Annual Review of Earth and Planetary Sciences, 48, 6.1-6.25.

535

https://doi.org/10.1146/annurev-earth-071719-054845

Confidential manuscript submitted to Journal of Geophysical Research: Solid Earth536

Saito, T., Baba, T., Inazu, D., Takemura, S., & Fukuyama, E. (2019). Synthesizing sea surface

537

height change including seismic waves and tsunami using a dynamic rupture scenario of538

anticipated Nankai trough earthquakes. Tectonophysics, 769, 228166.

539

https://doi.org/10.1016/j.tecto.2019.228166

540

Satake, K. (2002). Tsunamis. In W. H. K. Lee, P. Gennings, C. Kisslinger, & H. Kanamori

541(Eds.), International Handbook of Earthquake and Engineering Seismology (pp. 437–

542

451). London: Academic Press.

543

Scholz, C. H., & Campos, J. (1995). On the mechanism of seismic decoupling and back arc

544

spreading at subduction zones. Journal of Geophysical Research, 100(B11), 22103–

545

22115. https://doi.org/10.1029/95JB01869

546

Takahashi, N., Kodaira, S., Sato, T., Yamashita, M., Kaiho, Y., Miura, S., et al. (2015). Velocity

547

Structure of the Izu–Ogasawara Island Arc. Journal of Geography (Chigaku Zasshi),

548

124(5), 813–827 (in Japanese with English abstract).

549

https://doi.org/10.5026/jgeography.124.813

550

Titov, V., Rabinovich, A. B., Mofjeld, H. O., Thomson, R. E., & Gonzalez, F. I. (2005). The

551

global reach of the 26 December 2004 Sumatra tsunami. Science, 309, 2045–2048.

552

https://doi.org/10.1126/science.1114576

553

Tanioka, Y., & Satake, K. (1996). Fault parameters of the 1896 Sanriku tsunami earthquake

554

estimated from tsunami numerical modeling. Geophysical Research Letters, 23, 1549–

555

1552. https://doi.org/10.1029/96GL01479

556

Tsushima, H., Hino, R., Tanioka, Y., Imamura, F., & Fujimoto, H. (2012). Tsunami waveform

557

inversion incorporating permanent seafloor deformation and its application to tsunami

558

forecasting. Journal of Geophysical Research, 117, B03311.

559

https://doi.org/10.1029/2011JB008877

560
561

Uyeda, S., & Kanamori, H. (1979). Back-arc opening and the mode of subduction. Journal ofGeophysical Research, 84(B3), 1049–1061. https://doi.org/10.1029/JB084iB03p01049

562

Wallace, L. M., Webb, S. C., Ito, Y., Mochizuki, K., Hino, R., Henrys, S., et al. (2016). Slow slip

563

near the trench at the Hikurangi subduction zone, New Zealand. Science, 352, 701–704.

564

https://doi.org/10.1126/science.aaf2349

565
566

Webb, S. C. (1998). Broadband seismology and noise under the ocean. Reviews of Geophysics,
36, 105–142. https://doi.org/10.1029/97RG02287

Confidential manuscript submitted to Journal of Geophysical Research: Solid Earth567

Webb, S. C., & Nooner, S. L. (2016). High-resolution seafloor absolute pressure gauge568

measurements using a better counting method. Journal of Atmospheric and Oceanic

569

Technology, 33, 1859–1874. https://doi.org/10.1175/JTECH-D-15-0114.1

570

Wessel, P., Luis, J. F., Uieda, L., Scharroo, R., Wobbe, F., Smith, W. H. F., & Tian, D. (2019).

571

The generic mapping tools version 6. Geochemistry, Geophysics, Geosystems, 20,572

5556–5564. https://doi.org/ 10.1029/2019GC008515

573

Venkataraman, A., & Kanamori, H. (2004). Effect of directivity on estimates of a radiated574

seismic energy. Journal of Geophysical Research: Solid Earth, 109, B04301.

575

https://doi.org/10.1029/2003JB002548

576

Wilson, A., & Ma, S. (2021). Wedge plasticity and fully coupled simulations of dynamic rupture

577

and tsunami in the Cascadia Subduction zone. Journal of Geophysical Research: Solid

578

Earth, 126, e2020JB021627. https://doi.org/10.1029/2020JB021627

579

Ye, L., Lay, T., Kanamori, H., & Rivera, L. (2016). Rupture characteristics of major and great580(Mw ≥ 7.0) megathrust earthquakes from 1990 to 2015: 1. Source parameter scaling581

relationships. Journal of Geophysical Research: Solid Earth, 121, 826–844.

582

https://doi.org/10.1002/2015JB012426

Journal of Geophysical Research: Solid Earth
Supporting Information for
Earthquake rupture and tsunami generation of the 2015 Mw 5.9 Bonin event
revealed by in-situ pressure gauge array observation and integrated seismic andtsunami wave simulation
Tatsuya Kubota1, Tatsuhiko Saito1, Yoshio Fukao2, Hiroko Sugioka3, Aki Ito2, Takashi Tonegawa2,
Hajime Shiobara4, and Mikiya Yamashita5
National Research Institute for Earth Science and Disaster Resilience

1

Japan Agency for Marine-Earth Science and Technology

2

Kobe University

3

Earthqueke Research Institute, The University of Tokyo

4

Geological Survey of Japan, National Institute of Advanced Industrial Science and Technology

5Contents of this file
Texts S1 to S4
Figures S1 to S9Table S1Introduction
The procedure to synthesize the pressure waveform is described in Text S1. Text S2 presents aprocedure of the inversion for the coseismic slip distribution. The procedure and result for the
jackknife inversion test are described in Text S3. The procedure to calculate the radiatedenergy is described in Text S4. Figure S1 shows the data processing of the APG waveforms.
Comparison of the forward simulated hydrostatic pressure with and without considering thedispersion effect is shown in Figure S2. Figure S3 is a trade-off curve between the smoothing
weight and the RMSE value used to determine the smoothing weight in the inversion analysis.
The result for the jackknife inversion test is shown in Figures S4 and S5. The inversion results
using the Green’s function incorporating only the tsunami and only the dynamic pressurechange, are shown in Figures S6 and S7, respectively. Figure S8 depicts a comparison between
the observed waveform and the forward-calculated waveforms. Figure S9 shows the structure
model used for the forward calculation of the ocean-acoustic waves. Table S1 lists the locationof the APG stations used in this study.

1

Text S1.
This text explains the procedure to synthesize the APG waveforms containing staticand dynamic pressure changes. Although there is an approach to synthesize the pressurechange inside the focal area based on the fully coupled simulation of the seismic waves andtsunamis (e.g., Kozdon & Dunham, 2013; Lotto & Dunham, 2015; Madden et al., 2021; Maeda &Furumura, 2013; Saito et al., 2019), we adopted a simpler approach integrating conventional
calculation methods.
In this model, we assume pure-reverse faulting slip (rake = 90°) with a strike of 180°
for simplicity. We also assume a dip of 16° and the fault center depth measured from sea
surface at 16 km, which are determined based on the plate boundary model of Iwasaki et al.
(2015) and Takahashi et al. (2015). After some trial and error, we set the fault dimension of L =15 km and W = 15 km, and slip amount of D = 13.3 cm (seismic moment: M0 = 9.0 ×1017 Nm, μ= 30 GPa). We first calculate the static pressure change from the rectangular fault model. The
seafloor vertical displacement is calculated with an elastic half-space (Okada, 1992). In thiscalculation, we assumed the depth of the center of the rectangular fault is assumed at 16 km,
but the depth from the free surface (i.e., seafloor) of the half-space was set as 10 km
considering a water depth of 6 km (average depth of this region). We then incorporate thespatial filtering effect due to seawater of 6 km in depth for the seafloor displacement (Kajiura
1963; Saito, 2019) to calculate the sea-surface displacement. Using the sea-surface
displacement as an initial condition, we simulate a tsunami by solving the following linear
long wave equation (e.g., Satake, 2002):
!"
!#
!'
!#
!%
!#

!%

+ "$ ℎ$ !& = 0!%+ "$ ℎ$ !( = 0,
!"

(S1)

!'

= − !& − !(where the variable η is the sea surface height anomaly (tsunami height), M and N are thevelocity components integrated along the vertical direction over the seawater depth, h0 is thewater depth, and g0 =9.8 m/s2 is the gravitational constant. In this calculation, we used thebathymetry data published by the Japan Oceanographic Data Center (JODC), Japan Coast
Guard (JCG), with a spatial interval of Δx = 500 m. The temporal interval for the simulation was
set as Δt = 1 s. The sea-surface height change and seafloor vertical displacement are thenconverted to the static pressure change using the following relationship:')(*+,-#.#/0 ()) ≈ ,$ "$ -.()) − /1 ())0,(S2)

where ρ0 is the seawater density, and η(t) and uz(t) are the time history of the sea-surface
height (tsunamis) and seafloor vertical displacement, respectively. Here, we assumed that a 1
cm change of seawater column height (.()) − /1 ())) is converted to a hydrostatic pressurechange of 1hPa (i.e., ρ0g0 = 1 hPa/cm).

2

We then synthesize the dynamic pressure change. Assume a point source at thecenter of the rectangular fault with the same focal mechanism and the same seismic moment,
we calculate seafloor vertical accelerogram az(t) using a conventional elastodynamic equation
with discrete wavenumber method (e.g., Herrmann, 2013). To keep the consistency with theseafloor displacement calculation in the tsunami simulation, we use a simplified half-spacestructure model with VP = 7.00 km/s, VS = 4.04 km/s, and ρ = 2.45 g/cm3. The seafloor verticalaccelerogram is then converted to the dynamic pressure change using the following
relationship:'*(2.3/0 ()) ≈ ,$ ℎ$ 11 ()).

(S3)

Here, we assumed a seawater density of 1.03 g/cm3. The station depths (h0 in equation (S3))
are listed in Table S1 as well as in Fukao et al. (2021). We finally combine the synthesized staticand dynamic pressure waveforms:'()) = ')(*+,-#.#/0 ()) + '*(2.3/0 ()).

(S4)

We used the impulse response function in the calculations for the static anddynamic pressure changes described above. To incorporate the effects of the rupture
duration, the time function 2()) is convolved, as follows:3()) = ∫ 3/3456-7 () 8 )2() − ) 8 )5) 8 ,
9

;<#

!

!

2()) = : 61 − cos ; : <=,

(S5)
(S6)

where 3/3456-7 ()) is the synthesized waveform with the impulse response and Tr is the lengthof the base of the time function (i.e. rupture duration). We adopted Tr = 5 s. We finally apply
the same lowpass filter as applied to the observation.

3

Text S2.
In the inversion analysis, we adopt the following observation equation:
=>?
B
(S7)
>? A = ; < E,
CD@
where xobs is the data vector, G is the matrix consisting of the Green’s function, S is theconstraint matrix for spatial smoothing, and u is a vector representing the fault slip amount.
In the inversion, we assumed a planer fault with a dimension of 40 km length × 28
km width and divided it into 4 km × 4 km subfaults. The pressure changes from each sub-fault(Green’s function, G) are calculated using the same procedure as the forward simulation. In the
inversion, considering the filter cutoff (0.033 Hz) is sufficiently low for the M ~6 earthquake toneglect the temporal evolution of the rupture propagation across the fault, we assume theruptures of all sub-faults begin simultaneously (i.e., infinite rupture propagation velocity).
We use the weighting of the smoothing α of 0.3, which is determined based on the
trade-off between the smoothing weight and the root mean square error (RMSE, Figure S2).
The RMSE is defined as:
9

;

?@A
=>?
FGHI = J ∑'− L/ 0 ,/B9-L/
'

(S8)

-(2

where L/,C- and L/ is the i-th data of the observed and synthetic waveforms, and N is thetotal number of data samples. We also impose a nonnegativity constraint (Lawson & Hanson,1974).

4

Text S3.
This text explains the procedure and results to evaluate the uncertainty of theinversion analysis based on the jackknife, or leave-one-out, approach. Out of the eight APG
stations, we exclude one and use the other seven APG stations to estimate the mainshock slip
distribution (Figures S4a–S4h). The inversion setting is identical to the original one which used
all the APG stations. We then calculate the mean value (μ) and the standard deviation (σ) using
these inversion results (Figures S4i-S4j). The possible minimum and maximum slip amounts
are evaluated based on the slip distribution with μ−1σ and μ+1σ (Figures S4k–S4l).
As a result of the jackknife test, the constraint of the slip at the northeastern region
was not so good, when the B10 station was excluded (Figure S4h), which was located justabove the slip area and at the eastern edge of the APG array. This leads to a relatively highstandard deviation (Figure S4j). Meanwhile, when the mean slip and standard deviation were
calculated based on the inversions using the B10 data (Figures S4a–S4g), the mean slip wasvery similar to the original inversion result and the deviation was small (Figures S4m–S4p). Thismeans station B10 is important to constrain the slip distribution. Using the mean slip andstandard deviation with the jackknife inversions including the B10 data (Figures S4mS4p), theuncertainties of the seismic moment and the area of the slip with large slips (D > 0.2Dmax) were
estimated as 8.0–11.6 × 1017 Nm and 208–224 km2.
We further evaluate the uncertainty of the stress drop value, based on the results of
the jackknife test (Figure S5). The stress drop distribution calculated by each slip distributionof the jackknife inversions are shown in Figures S5a–S5h. We then calculate the mean stressdrop μ and standard deviation σ, using all the jackknife inversion models (Figures S5i–S5l) andusing the jackknife inversion models including the B10 data (Figures S5m–S5p). Based on the
stress drop distributions with μ−1σ and μ+1σ, calculated from the inversions using the B10
data (Figures S5o and S5p), the uncertainty of the stress drop is estimated as 0.44–0.50 MPa.

5

Text S4.
This text explains the procedure to estimate the radiation energy ER and radiation
efficiency ER/M0 (Venkataraman & Kanamori, 2004). The radiation energy ER is expressed as(Venkataraman & Kanamori, 2004):
;

I
E<
E<
Ṙ (O)P 5O,
ID = M9FGH # + 9$GH #N ∫$ O ; PG
"

$

(S9)

Ṙ (O) is the moment rate
where ρ is the density, VP and VS are the P and S wave velocities, G
spectrum. We here assume the ω−2 model (Aki, 1967), for the shape of the moment rate
spectrum:'Ṙ (O) = "% J& ,
G
'
'
J KJ&

(S10)

where fc is the corner frequency. We set the corner frequency based on the optimum duration
estimated from the forward simulation of the ocean-acoustic waves, Tr = 7 s (Figure 4), as fc =1/Tr = 0.14 Hz. We use the density and P and S wave velocities from the values used for theforward simulation of the dynamic pressure (see Text S1). We used 0.01 to 1 Hz for the
integration. As a result, we obtain ER = 1.77 ×1012 Nm and ER/M0 = 1.87 ×10−6.
We also evaluate the uncertainty of the radiated energy and radiation efficiency.
Based on the numerical simulation of the ocean-acoustic waves (Figure 4), the simulations
with the duration Tr = 6 or 8 s reasonably explained the observed pressure waveforms as wellas Tr = 7 s. When assuming Tr = 6 s and 8 s, we obtain ER = 2.70 ×1012 Nm and 1.23 ×1012 Nm,
respectively. Using the smallest value of the possible range of the moment magnitude (M0 =8.0 × 1017 Nm, Text S3) and the radiation efficiency with Tr = 6 s, the upper limit of the radiation
efficiency was estimated as ER/M0 = 3.33 ×10−6. Even if taking account of these uncertainties,
these values are consistent with those for tsunami earthquakes estimated from the global M >~7 event catalog (Ye et al., 2016).

6

(a) Detided pressure(b) Processed pressure

B05.H1

LPF(100s)
0.6cm/s

B05.BZ

BPF(30-100s) LPF(30s)

B01
B02
B03
B06
B07

B05.H2

B05.DPG

B01

B08

B02
B03
0.1MPa ∆P

B06

B09

B07
5hPa ∆P

B08
B10

B09
B10
−600

Foreshock
Mainshock

Foreshock
Mainshock06001200 −120 −60
0
60
120
180
240
300
Time [s]Time [s]Figure S1. Data processing of the APG records. (a) De-tided waveforms. Seismograms recorded by BBOBS (instrument responses
are not corrected) and the pressure recorded by differential pressure gauge (DPG, response-corrected) at station B05 are also shown.
(b) Filtered waveforms. Blue dashed, green dashed-and-dotted, and red solid traces denote the lowpass-filtered waveforms (a cutoff
of 0.01 Hz), bandpass-filtered waveforms (0.01 to 0.033 Hz), and the lowpass-filtered waveforms (a cutoff of 0.033 Hz), respectively.
Note that t = 0 s in the horizontal axis is the origin time. The focal times of the mainshock and foreshocks are also marked. Note thatthe filters were applied in both the forward and backward directions.
1Observed
Without dispersion

With dispersion

B01
B02
B03
B06
B07
B08

B09

5hPa ∆P

B10

Lowpass 0.033 Hz
−120

−60

0

60
120
Time [s]

180

240

300

Figure S2. Comparison of the observed (black) and simulated hydrostatic pressure
waveforms from the rectangular fault in Figure 1a, considering the dispersion effect ornot. Blue traces are the synthetic hydrostatic pressure waveforms without thedispersion effect (identical to those in Figure 1d), and green traces are the waveforms
including the dispersion effect. The simulated waveforms are very similar to eachother, indicating the effect of the dispersion is very minor.

1

Figure S3. Trade-off curve between the smoothing weight α and the root mean squareerror (RMSE) for the inversion considering only static (blue), only dynamic (green), andboth dynamic and static pressure changes (red). In the inversion, we adopted α = 0.3.

2

141.5°142°142.5°(a) Without B01

141.5°

142°

142.5°(b) Without B02

141.5°

142°

142.5°(c) Without B03

141.5°

142°

142.5°(d) Without B06

31.5°

B07 B08 B09 B10B07 B08 B09 B10B06

B07 B08 B09 B10B06

B02 B03

B02

B03

31°

B01
0

10

20

M0= 9.5e+17Nm

[cm]
0

5

10

15

20

0

B02 B03

B01
10

20

M0= 9.8e+17Nm(e) Without B07

B07 B08 B09 B10

B06

[cm]
0

5

10

15

20

0

B01
10

20

M0= 9.5e+17Nm

[cm]
0

5

10

15

20

(g) Without B09(f) Without B08

0

10

20

10

20

10

20

M0= 1.0e+18Nm

[cm]
0

5

10

15

20

0

5

10

15

20

0

5

10

15

20

10

15

20

(h) Without B10

31.5°

B08 B09 B1031°

B09 B10

B07

B06

B07 B08

B06

B06

B02 B03B02 B03B02 B03

B01

B01

B01

0

10

20

M0= 1.0e+18Nm

[cm]
0

5

10

15

20

0

10

20

M0= 1.1e+18Nm(i) Mean µ (from (a)-(h) models)

B07 B08 B09

B10

B06

[cm]
0

5

10

15

20

0

B02 B03B01
10

20

M0= 8.2e+17Nm

[cm]
0

5

10

15

20

(k) µ - 1σ(j) Standard deviation σ0

M0= 1.1e+18Nm[cm](l) µ + 1σ

31.5°

31°
0

10

20

M0= 9.9e+17Nm

[cm]
0

5

10

15

20

(m) Mean µ (from (a)-(g) models, with B10)

0

10

20

[cm]
0

5

10

15

20

0

10

20

M0= 5.5e+17Nm(n) Standard deviation σ (from (a)-(g) models)

[cm]
0

5

10

15

20

0

M0= 1.4e+18Nm(o) µ - 1σ (from (a)-(g) models)[cm](p) µ + 1σ (from (a)-(g) models)

31.5°

31°
0

M0= 9.8e+17Nm

10

20

[cm]
0

5

10

15

20

0

10

20

[cm]
0

5

10

15

20

0

M0= 8.0e+17Nm

10

20

[cm]
0

5

10

15

20

0

M0= 1.2e+18Nm

10

20

[cm]
0

5

Figure S4. Result of the jackknife inversion test. (a–h) Slip distributions based on the inversions without stations (a) B01, (b) B02,(c) B03, (d) B06, (e) B07, (f) B08, (g) B09, and (h) B10. (i–l) Mean slip μ, standard deviation σ, distribution of μ+1σ and μ−1σ, basedon all jackknife inversion models ((a)–(h) models), respectively. (m–p) Same as Figure S4i–S4l, but based on the jackknife inversion
models using B10 ((a)–(g) models).
1141.5°142°142.5°(a) Without B01

141.5°

142°

142.5°(b) Without B02

141.5°

142°

142.5°(c) Without B03

141.5°

142°

142.5°(d) Without B06

31.5°

B07 B08 B09 B10B07 B08 B09 B10B06

B07 B08 B09 B10B06

B02 B03

B02

B03

31°

B01
0

10

∆σE = 0.51MPa

20−1.0 −0.5 0.0 0.5 1.0

[MPa]

0

B02 B03B01
10∆σE = 0.49MPa(e) Without B07

B07 B08 B09 B10B06

20−1.0 −0.5 0.0 0.5 1.0

[MPa]

0

B01
10

∆σE = 0.50MPa20−1.0 −0.5 0.0 0.5 1.0[MPa](g) Without B09(f) Without B08

0

10

20
−1.0 −0.5 0.0 0.5 1.0

0

10

20
−1.0 −0.5 0.0 0.5 1.0

0

10

20
−1.0 −0.5 0.0 0.5 1.0∆σE = 0.50MPa[MPa](h) Without B10

31.5°

B08 B09 B1031°

B09 B10

B07

B06

B07 B08

B06

B06

B02 B03B02 B03B02 B03

B01

B01

B01

0

10

∆σE = 0.48MPa

20−1.0 −0.5 0.0 0.5 1.0

[MPa]

0

10

∆σE = 0.46MPa(i) Mean µ (from (a)-(h) models)

B07 B08 B09

B10

B06

20
−1.0 −0.5 0.0 0.5 1.0

[MPa]

0

B02 B03B01
10∆σE = 0.37MPa

20−1.0 −0.5 0.0 0.5 1.0[MPa]∆σE = 0.33MPa(k) µ - 1σ(j) Standard deviation σ[MPa](l) µ + 1σ

31.5°

31°
0

10

∆σE = 0.41MPa

20−1.0 −0.5 0.0 0.5 1.0[MPa](m) Mean µ (from (a)-(g) models, with B10)

0

10

20
−1.0 −0.5 0.0 0.5 1.0

[MPa]

0

10

∆σE = 0.39MPa(n) Standard deviation σ (from (a)-(g) models)20−1.0 −0.5 0.0 0.5 1.0[MPa]∆σE = 0.41MPa(o) µ - 1σ (from (a)-(g) models)[MPa](p) µ + 1σ (from (a)-(g) models)

31.5°

31°
0

∆σE = 0.47MPa

1020−1.0 −0.5 0.0 0.5 1.0

[MPa]

0

10

20
−1.0 −0.5 0.0 0.5 1.0

[MPa]

0

∆σE = 0.46MPa

1020−1.0 −0.5 0.0 0.5 1.0

[MPa]

0

∆σE = 0.50MPa

1020−1.0 −0.5 0.0 0.5 1.0[MPa]

Figure S5. Result of the jackknife inversion test. (a–h) Stress drop distributions based on the inversions without stations (a) B01,(b) B02, (c) B03, (d) B06, (e) B07, (f) B08, (g) B09, and (h) B10. (i–l) Mean stress drop μ, standard deviation σ, distribution of μ+1σ and
μ−1σ, based on all jackknife inversion models ((a)–(h) models), respectively. (m–p) Same as Figures S5i–S5l, but based on the
jackknife inversion models using B10 ((a)–(g) models)2141.5°142°142.5°

141.5°

142°

142.5°(b)Stress drop(a)Slip distribution
31.5°

31.5°
Model

B07 B08 B09 B10B07 B08 B09 B10

B06

B06

B02 B03B02 B03
31°

0

Pressure Change [hPa]

40

B01

31°

B01
10

20

[cm]

Contour:1cm

0

5

10

15

(c) Inverted

Synthesized (Static Pressure)020(d) Observed

B01
B02
35 B03
B06
B07
30
B08

10

Dynamic20−1.0 −0.5 0.0 0.5 1.0[MPa]

Static+Dynamic

B01
B02
B03
B06
B07
B08

25
B09

B09

20
15
10

5hPa ∆P

B10

B10

5
0
−120

−60

0

60
120
Time [s]180240300 −120

−60

0

60
120
Time [s]

180

240

300

Figure S6. Result of the fault slip inversion using the Green’s function incorporating
only static pressure change. See Figure 3 for the caption of this figure.

1

141.5°142°142.5°

141.5°

142°

142.5°(b)Stress drop(a)Slip distribution
31.5°

31.5°
Model

B07 B08 B09 B10B07 B08 B09 B10

B06

B06

B02 B03B02 B03
31°

0

Contour:1cm

Pressure Change [hPa]

40

B01

31°

B01

(c) Inverted

10

20

[cm]
0

5

10

15

0

20

Synthesized (Dynamic Pressure)

B01
B02
35 B03
B06
B07
30
B08

(d) Observed

Static1020−1.0 −0.5 0.0 0.5 1.0[MPa]

Static+Dynamic

B01
B02
B03
B06
B07
B08

25
B09

B09

20
15
10

5hPa ∆P

B10

B10

5
0
−120

−60

0

60
120
Time [s]180240300 −120

−60

0

60
120
Time [s]

180

240

300

Figure S7. Result of the fault slip inversion using the Green’s function incorporating
only dynamic pressure changes. See Figure 3 for the caption of this figure.

2

Pressure Change [hPa]

40

Observed
Mainshock
Foreshock (0.1×Mainshock) Mainshock+Foreshock

B01
B02
35 B03
B06
B07
30
B08
25
B09
20
15
10

5hPa ∆P

B10

5
0
−120

−60

0

60
120
Time [s]

180

240

300

Figure S8. Comparison between the observed waveform (gray) and the forwardcalculated waveforms for foreshock (green), mainshock (blue), and both foreshock andmainshock (red), based on the mainshock slip distribution (Figure 3). The red lines inthis figure is identical to those in fig. 4 of Fukao et al. (2021).

3

Figure S9. Velocity structure for simulation of the ocean-acoustic wave.

4

Table S1. Location of the array observation instrumentsa
StationInstrumentb
Longitude [°E]
Longitude [°N]
Depth [m]B01
APG
141.7118
30.9985
5313
B02
APG
141.6590
31.0759
5155
B03
APG
141.7630
31.0754
5601
B04
APG c
141.819431.1530
4794
B05
BBOBS+DPG
141.7178
31.1536
5427
B06
APG
141.6056
31.1541
4853
B07
APG
141.5536
31.2322
4671
B08
APG
141.6594
31.2319
5089
B09
APG
141.7694
31.2343
5077
B10
APG
141.8699
31.2298
5670
a
See Fukao et al. (2021) for more details of the observation. Sampling rate of the original APG
data is 4 Hz, but the data are sampled to 1 Hz for the present analyses.
b
APG: Absolute pressure gauges, DPG: Differential pressure gauges, BBOBS: Broadband oceanbottom seismometer.
c
No data was obtained due to an error of data logging.

5