manuscript submitted to Journal of Geophysical Research: Solid Earth1

Improving the constraint on the Mw 7.1 2016 off-Fukushima shallow normal-

2

faulting earthquake with the high azimuthal coverage tsunami data from the3

S-net wide and dense network: Implication for the stress regime in the4

Tohoku overriding plate5
6

Tatsuya Kubota1, Hisahiko Kubo1, Keisuke Yoshida2, Naotaka Y. Chikasada1, Wataru

7

Suzuki1, Takeshi Nakamura3, and Hiroaki Tsushima4891National Research Institute for Earth Science and Disaster Resilience, Tsukuba, Japan.

102Graduate School of Science, Tohoku University, Sendai, Japan.

113Central Research Institute of Electric Power Industry, Abiko, Japan.

124Meteorological Research Institute, Japan Meteorological Agency, Tsukuba, Japan

13
14

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

15
16
17

Key Points:•

observed by the S-net wide and dense pressure gauge network

18
19

•

22
23

Use of the near-field and the high-coverage array significantly improved the constraint of
the fault modeling of the 2016 earthquake

20
21

Tsunamis due to the 2016 off-Fukushima shallow normal-faulting earthquake were•

Horizontal extensional stress predominant even before the 2011 Tohoku earthquake
should be the main cause of the earthquake and tsunami

manuscript submitted to Journal of Geophysical Research: Solid Earth

24
25

Abstract
Tsunamis with amplitudes of up to 40 cm, related to the Mw 7.1 normal-faulting

26

earthquake off Fukushima, Japan, on November 21, 2016, were clearly recorded by a new27

offshore wide and dense ocean-bottom pressure gauge network, S-net, with high azimuthal

28

coverage located closer to the focal area. We processed the S-net data and found some stations

29

included the tsunami-irrelevant drift and step signals. We analyzed the S-net data to infer the30

tsunami source distribution. A subsidence region with a narrow spatial extent (~40 km) and a31large peak (~200 cm) was obtained. The other near-coastal waveforms not used for the inversion

32

analysis were also reproduced very well. Our fault model suggests the maximum stress drop

33

across the fault plane of > ~10 MPa and the average of 4.2 MPa, whereas the shear stress34

increase along the fault caused by the 2011 Tohoku earthquake was only ~2 MPa. Past studies

35

have suggested that horizontal compressional stress around this region switched to horizontal

36

extensional stress after the Tohoku earthquake due to its stress perturbation. The present result,37

however, suggests that the horizontal extensional stress was locally predominant at the38

shallowest surface around the focal area even before 2011. The present study demonstrates that39

the S-net high-azimuthal-coverage pressure data provides a significant constraint on the fault

40

modeling, which enables us to discuss the stress regime within the overriding plate at the41

offshore. Our analysis provides an implication for crustal stress states, which is important for42

understanding generation mechanisms of intraplate earthquakes.

43
44
45

Plain Language Summary
On November 21, 2016 (UTC), a large earthquake occurred within the continental plate46

off Fukushima, Japan, and a new seafloor tsunami network, S-net, recorded its tsunamis with47

much higher azimuthal coverage and with shorter epicentral distance than any of the previous

48

networks. We analyzed the S-net data to reveal the rupture process of this earthquake. Our result

49

explained all of the S-net data and the other tsunami network data very well. According to past50

studies, the continental plate in northeastern Japan was under horizontal compression before the51

2011 Tohoku earthquake due to the pushing force by the subducting oceanic plate. However, our

52

rupture modeling result suggested that the plate around the earthquake rupture area was53

horizontally stretched even before the Tohoku earthquake, so that the off-Fukushima earthquake54

occurred. Our study demonstrated that the S-net, which has high spatial coverage, makes itmanuscript submitted to Journal of Geophysical Research: Solid Earth55

possible to reveal the rupture model of offshore earthquakes, which was difficult in the past

56

before S-net became available. The S-net will also enable us to discuss the impact of the Tohoku

57

earthquake on the crustal stress, which is necessary for understanding the earthquake generation

58

mechanics.

59

manuscript submitted to Journal of Geophysical Research: Solid Earth60

1 Introduction

61

In this decade, the coseismic rupture process of the 2011 Tohoku earthquake and its

62

preseismic and postseismic processes have been investigated in detail (e.g., Hino, 2015; Kodaira63

et al., 2020; 2021; Lay, 2018; Uchida & Bürgmann, 2021; Wang et al., 2018). In response to the64

Tohoku earthquake, a new wide offshore deep-ocean observation network, Seafloor Observation

65

Network for Earthquakes and Tsunamis along the Japan Trench (S-net), has been constructed off

66

eastern Japan (Aoi et al., 2020; Kanazawa et al., 2016; Mochizuki et al., 2017; Uehira et al.,67

2016, Figure 1a). Recent studies have started to utilize S-net ocean-bottom seismometers to68

investigate the seismotectonics and geodynamics in the Tohoku subduction zone (Dhakal et al.,692021; Hua et al., 2020; Matsubara et al., 2019; Nishikawa et al., 2019; Sawazaki & Nakamura,

70

2020; Takagi et al., 2019, 2021; Tanaka et al., 2019; Uchida et al., 2020; Yu & Zhao, 2020). The

71

S-net also incorporates ocean-bottom pressure gauges (OBPGs), which are expected to be72

utilized for tsunami forecasts (e.g., Aoi et al., 2019; Inoue et al., 2019; Mulia & Satake, 2021;

73

Tanioka, 2020; Tsushima & Yamamoto, 2020; Wang et al., 2021; Yamamoto et al., 2016a;

74

2016b). The other potential contributions to the earth sciences of the S-net OBPG have also been75

demonstrated, such as understanding the wave propagation process in the ocean as well as the76

rupture process of subseafloor earthquakes (Kubota et al., 2020a; 2021; Saito & Kubota, 2020;

77

Saito et al., 2021). The wide and dense network data of S-net will significantly broaden our

78

understanding of the Tohoku subduction zone after the Tohoku earthquake.

79

On November 21, 2016, a major shallow normal-faulting earthquake occurred within

80

the overriding plate off Fukushima Prefecture (20:59 UTC, Mw 6.9, 12 km, Global CMT

81[GCMT], https://www.globalcmt.org, Figure 1, hereafter referred to as the off-Fukushima

82

earthquake). Compared with the GCMT centroid, its epicenter, as determined by Japan

83

Meteorological Agency (JMA), was located ~20 km east to northeast (white star in Figure 1).

84

Numerous aftershocks accompanied this earthquake (Figures 1b and 1c). It has been reported

85

that the tsunamis associated with the off-Fukushima earthquake were observed by onshore and86

offshore tsunami networks (e.g., Gusman et al., 2017; Kawaguchi et al., 2017; Suppasri et al.,872017). However, these stations were located only on the shore-side from the focal area, and the88source-station distances are large (Figure 1a). In contrast, the S-net OBPGs recorded tsunamis

89

with much higher azimuthal coverage and with a closer distance to the focal area (~30 km,

90

Figure 1a). Because of the much better station coverage of the S-net, the constraint on the initialmanuscript submitted to Journal of Geophysical Research: Solid Earth91

sea height (tsunami source) estimation and the finite fault modeling of the off-Fukushima

92

earthquake will be significantly increased, as compared with the previous datasets.

93

The normal-faulting mechanism of the off-Fukushima earthquake is similar to nearby

94

shallow normal-faulting micro-seismicity within the overriding plate, with a tensile axis (σ3)

95oriented in basically the east-west direction, which significantly increased after the Tohoku

96

earthquake (Figures 1d–1f, e.g., Asano et al., 2011; Hardebeck & Okada, 2018; Hasegawa et al.,97

2012; Tanaka et al., 2014; Wang et al., 2019; Yoshida et al., 2012). This increase in the normal-

98

faulting seismicity is considered to be related to the significant stress perturbation by the Tohoku

99

earthquake, which switched the intraplate stress regime from horizontal compression to100

horizontal extension (e.g., Hasegawa et al. 2012). If we can obtain a detailed fault model of the101

off-Fukushima earthquake, then the quantitative relationship between the crustal stress released

102

during the off-Fukushima earthquake (i.e., stress drop) and the stress increase due to the 2011

103

Tohoku earthquake can be discussed. Because the generation of earthquakes is very closely104

related to the process of the stress accumulation and release, this quantitative comparison of the105

stresses is essential to deepen our understanding of the temporal change of the crustal stress state106

associated with the Tohoku earthquake and our knowledge of the generation mechanisms of the107

subseafloor crustal earthquakes related to the offshore megathrust earthquake, which may excite

108

the significant tsunamis to cause a severe damage to the coast.

109

Because of the tsunamis' much smaller propagation velocity than that of seismic110

waves, the tradeoff between the earthquake source dimension and the rupture velocity is much

111

smaller (Kubota et al. 2018a). In addition, shallow earthquakes generally excite tsunamis more112

efficiently than deep earthquakes (Kubota et al. 2019). So tsunami data have a strong advantage

113

in the robust constraint on the fault slip extent and slip amount of the 2016 off-Fukushima

114

earthquake, leading to the robust constraint on the stress drop. In the present study, therefore, we

115

estimate the detailed finite fault model of the off-Fukushima earthquake using the S-net OBPG

116

data. From the finite fault model, we also attempt to examine the normal-faulting stress state

117

within the crust around the off-Fukushima earthquake and its relationship with the Tohoku

118

earthquake, based on the stress drop estimation from the finite fault model. Section 2 describes

119

the dataset used in this study, and Section 3 summarizes the feature in the S-net OBPG data. The

120

spatial distribution of the initial sea surface height (tsunami source) and the finite fault model of121

the off-Fukushima earthquake are estimated in Sections 4 and 5, respectively. Section 6

manuscript submitted to Journal of Geophysical Research: Solid Earth122

examines the relationship between the Tohoku earthquake and the stress regime around the focal

123

area. Section 7 concludes the present study.

124(b) 2016.11.21EQ - 2016.11.30

aki
Ibar
207

35°
140°

40

0 20 40140.5°

141°141.5°

142°

20

Depth [km]

50

30

0
10
20
30
40
50

A’

0

AoA70

10

20 G19

20

21
26 25
24 23 22
01 02
03
04
08
05
09 07 06
11 10
12

S2
S1

141°

142°

37°

36.5°
143°

142.5°FKOB2
16 17 18 FKOB1
19

209

37.5°

(c)

AoA60(G17)

AoA50

206 15

36°

G01

38°

0

100

4000

8000

0

Sea depth [m]
143°

144°

145°

30
40

A(NW)
−40

A’(SE)−20

020Depth [km]

37°

A

60

Iw

38°

80

S4
S3

70

23
28 27 26 25
24
01
02 03
04
05
TM2
TM1
06
11 10
12
07
08
09
803
13
14 15
218
16 17
801
18
G12 19205 26 25 24 23
22 21 2001
214
02 03 04
05
G15
06
07
Fukushima
G16
13 12
08
11
10 09
806 14

Mi
ya
gi

39°

Depth < 50 km

(a)

ate

40°

50

40

Horizontal Distance [km]38°(d) 1997-2011.3.11(e) 2011.3.11EQ-2016.11.21EQ (f) 2016.11.21EQ-2017.12.31
Depth < 20 km

37°

125

141°

142°

143°

141°

142°

143°

141°

142°

143°

126

Figure 1. (a) Location map of the present study. Locations of the tsunami stations are shown by

127

colored symbols (black circle: S-net OBPG, blue inverted triangle: ERI OBPG, yellow inverted

128

triangle: Tohoku University OBPG, orange square: NOWPHAS GPS buoy, pink triangle:

129

NOWPHAS wave gauge). The epicenter (white star) and the CMT solution (red) of the off-130

Fukushima earthquake are taken from JMA and GCMT, respectively. (b) Enlarged view of the131

rectangular area drawn by gray lines in Figure 1a. Aftershocks during about one week as

132

determined by JMA are shown (color denotes its depth). Orange contours show the depth of the133

subducting plate interface (Nakajima & Hasegawa, 2006). The locations of fresh seafloor cracks

manuscript submitted to Journal of Geophysical Research: Solid Earth

134

found by the JAMSTEC survey are shown by blue triangles. (c) Vertical cross section along line

135

A-A' in Figure 1b. (d–f) The F-net fault mechanisms (Fukuyama et al., 1998) at depths shallower

136

than 20 km, (d) before the Tohoku earthquake, (e) between the Tohoku earthquake and the off-137

Fukushima earthquake, and (f) after the off-Fukushima earthquake.

138
139

2 Tsunami dataset
The present study used the S-net OBPG data (black circles in Figure 1a, Wang &

140
141

Satake, 2021). Although S-net now consists of 150 observatories (Aoi et al., 2020), 25 of these142

observatories, located at the outer-trench region, were not installed when the off-Fukushima

143

earthquake occurred. Each observatory is equipped with absolute pressure sensors manufactured

144

by Paroscientific, Inc. (e.g., Polster et al., 2009; Watts & Kontoyiannis, 1996). Two pressure

145

sensors are equipped in each observatory for redundancy. The sensors are not directly exposed to146

the seawater, but rather are sealed in a metal housing filled with oil. The metal housing is further

147

sealed in a metal cylindrical vessel filled with oil. The external pressure is transferred to the148pressure sensor inside via a diaphragm made of hard rubber. See Aoi et al. (2020) for more

149

details.

150

In addition to S-net, we use other OBPGs to evaluate the modeling resolution. We use

151

the OBPGs off Iwate Prefecture installed by the Earthquake Research Institute (ERI) of the152

University of Tokyo (blue inverted triangles in Figure 1a, Gusman et al., 2017; Kanazawa &153

Hasegawa, 1997) and the OBPGs off eastern Japan installed by Tohoku University (yellow

154

inverted triangles, Hino et al., 2014; 2021). We also use the offshore GPS buoys (orange

155

squares) and wave gauges (pink triangles) of the Nationwide Ocean Wave information network

156

for Ports and HArbourS [NOWPHAS] (Kawaguchi et al., 2017; Nagai et al., 1998).

157
158
159

3 Fundamental feature of the S-net OBPGs: Tsunami-irrelevant pressure signals
In order to investigate the fundamental feature of the S-net OBPG signals, we first160

process the OBPG data. We decimate the original 10 Hz data to 1 Hz (Figures 2a–2i). We then161

subtract the theoretical tide calculated by the model of Matsumoto et al. (2000) and apply a low-162

pass filter with a cutoff of 100 s in both forward and backward directions, to reduce the high-163

frequency seismic wave signals (Figure 2j–2r).

manuscript submitted to Journal of Geophysical Research: Solid Earth164

Figures 2a–2i show the 1-Hz-sampling pressure waveforms. The high-frequency

165

fluctuations related to the seismic waves and ocean-acoustic waves (e.g., Kubota et al., 2020b)

166

are observed. The gradual pressure increases related to the ocean tide are also observed, although167

some traces show different trends. The pressure changes recorded by the two sensors equipped in168

the same observatory (black and gray lines) are very similar to each other. The difference

169

between these two traces (red lines) is around zero, although some stations have offsets in the170

differences. At station S2N13, which is located just above the focal area of the off-Fukushima

171

earthquake, no seismic or tsunamis signals were recorded, although the co-equipped seismometer

172

correctly recorded the ground shaking of this earthquake (Dhakal et al., 2021; Takagi et al.,173

2019; see also https://www.hinet.bosai.go.jp/topics/off-fukushima161122/?m=snet, in Japanese).

174

This may suggest that the pressure observation part at S2N13 observation node did not work

175

correctly. We note that the instrument at this site was replaced in 2020176

(https://www.seafloor.bosai.go.jp/notice/notice_200414_1.pdf, in Japanese), and this site

177

recorded tsunamis related to a MJMA 7.3 earthquake on 13 February 2021

178

(https://www.hinet.bosai.go.jp/topics/off-fukushima210213/?m=others, in Japanese).

179

180

manuscript submitted to Journal of Geophysical Research: Solid Earth181

Figure 2. (a–i) The 1-Hz ocean-bottom pressure waveforms for stations (a) S2N02, (b) S2N03,

182(c) S2N06, (d) S2N11, (e) S2N12, (f) S2N13, (g) S2N14, (h) S2N15, and (i) S2N18. Black and183

gray traces denote the waveforms from each of the pressure sensors. Red traces denote the184

difference between the two sensors. Note that the vertical scale for the difference waveforms is185

different in each subfigure. The dominant signals are indicated by arrows and text. The epicentral

186

distance Δ measured from the JMA epicenter, and the tilt change λ, rotation angle change θ, and187

peak ground acceleration (PGA) values measured by the co-equipped accelerometer (Takagi et188al., 2019) are also shown. (j–r) Ocean-bottom pressure waveforms after data processing for the189stations.

190
191

Although the tsunamis are confirmed in the lowpass-filtered S-net OBPG waveforms

192(Figures 2j–2r), we also recognize some signals irrelevant to the tsunamis, such as the large drift

193

components (e.g., > ~50 hPa/hour, S2N06 and S2N18). If these drifts are caused by a vertical194

movement due to the postseismic deformation or the unstable slow sliding of the sensor, the195

vertical movement rate will be ~50 cm/hour (1 hPa pressure change is approximated as 1 cm

196

displacement). This is too unrealistic. In addition, these drift features are also confirmed in some197of the other sites at water depth shallower than 1500 m, in which the instruments are buried

198

about one meter beneath the seafloor (Aoi et al., 2020). Therefore, these drifts are not caused by199

the real movement, but it might be possible to be caused by the mechanical reason. Some

200

previous studies reported that the Paroscientific pressure sensors contain instrumental drift with

201

rates of ~8.8 hPa/year (Inazu & Hino, 2011; Polster et al. 2009; Watts & Kontoyiannis, 1996),

202

although this previously-reported rates are much smaller. It is also incomprehensible that the203

drift rates are completely identical in the sensor pair, although the instrumental drift must be

204

individually different in each sensor. Therefore, we do not consider the cause of these drifts to be205

the one previously reported. Although we cannot identify the reason for these drifts, we suspect206

the observation system of the S-net may be relevant. The observation system of the S-net

207

observation node, which includes not only the OBPG sensors but also other instruments such as

208

seismometers, a power supply unit, and a real-time data transmitting unit (Aoi et al., 2020), is

209

much more complicated than the ordinary offline pop-up pressure observation (e.g., Hino et al.,2102014).

manuscript submitted to Journal of Geophysical Research: Solid Earth211

In addition, abrupt steps at the origin time are observed at some OBPGs, particularly at

212

S2N11, S2N12, S2N14, and S2N15. The step is also observed at S2N13, where no tsunami

213

signals were recorded. If we consider the pressure offset changes as a result of the seafloor

214

vertical movement, these pressure changes correspond to a seafloor vertical displacement of

215~30–60 cm (1 hPa pressure change is approximated as 1 cm vertical movement). Considering the216

source-station distances, these displacements seem too large compared with those expected from217

typical M~7 earthquakes. Furthermore, even if the OBPGs are located inside the focal area

218

where the vertical displacement is large, the ocean-bottom pressure, or the seawater column

219

height above the OBPG, cannot change so abruptly because both seafloor and sea-surface

220

simultaneously move vertically during tsunami generation (e.g., Tsushima et al. 2012).

221

Therefore, these steps are unlikely to be caused by the seafloor permanent displacement. Similar222

pressure steps were also recorded by the S-net and the other OBPG networks during the past

223

earthquakes (Kubota et al., 2018b; 2020a; Wallace et al. 2016), which are not considered to be224

related to the tsunami or the seafloor crustal deformation.

225

It has been reported that outputs of Paroscientific pressure sensors strongly depend on

226

its orientation relative to the direction of gravity (Chadwick et al., 2006). Thus, the step signals

227

might be caused by the rotation of the pressure sensor. According to Chadwick et al. (2006), the228

rotation angle change of the pressure sensor of θ ~10° roughly corresponds to the apparent229

pressure offset change of up to ~10 hPa. Takagi et al. (2019) analyzed the co-equipped230

accelerometer during the off-Fukushima earthquake and found that some observatories near the231

epicenter rotated associated with large seafloor ground motion (Figure 2). However, comparing232

the rotation angles at some near-source stations (e.g., θ = 0.86° at S2N12 and 9.95° at S2N14,

233

Takagi et al., 2019), the observed pressure steps were extremely large (> ~50 hPa). Furthermore,234

considering that the sensitivity to the rotation angle must be different in each sensor, it is quite

235

strange that the amounts of the pressure step in two pressure sensors are almost identical. We

236

also confirm that the pressure steps in the two pressure sensors are different at some stations

237

where the large rotation was observed (e.g., S2N13, S2N15), leading to the steps around the238

focal time in the difference traces between the two sensor outputs (red lines in Figure 2). Taking

239

these points into account, we consider that the dominant cause of the pressure steps is not the240

response to the sensor rotation as reported by Chadwick et al. (2006), and the difference in the241

steps between the two sensors may be due to the difference in the response to the rotation angle.

manuscript submitted to Journal of Geophysical Research: Solid Earth242

Although more detailed investigation is needed, we suspect the strong shaking of the instrument

243

due to the seismic motion might have mainly caused these steps, which may affect the244

observation system of S-net, such as the transferring system of the external pressure to the sensor

245

inside of the metal housing. As a summary of this section, we emphasize that we must be careful

246

to analyze the OBPG data to distinguish whether such signals are real or are artifacts related to247

the drift or offset, although the S-net OBPGs clearly recorded the tsunamis due to the 2016 off-248

Fukushima earthquake.

249
250

4 Tsunami source modeling251

4.1 Modeling procedure252

In this section, we analyze the S-net data to estimate the spatial distribution of initial253

sea-surface height (tsunami source) of the off-Fukushima earthquake and to investigate how the254

S-net OBPGs provide better constraint. In order to reduce the long-period tsunami-irrelevant

255

drift signals as well as the short-period seismic wave components, we apply the bandpass filter

256

with passbands of 100–3,600 s (Figure 3b). We here briefly describe the procedure for the257

tsunami source modeling. The full details are shown in Text S1.

258

We distribute the unit source elements of the seafloor vertical displacement with259

horizontal spatial intervals of 2 km, in an area of 50 km × 50 km (rectangular area in Figure 3a).

260

To calculate pressure change waveforms excited by each unit source element (i.e., the Green's

261

functions) we simulate tsunamis by solving a linear dispersive tsunami equation (Saito, 2019;262

Saito et al., 2010). We use the JTOPO30 bathymetry data with a spatial resolution of 30 arcsec,

263

interpolating the spatial interval of 1 km. We assume the vertical displacement of the unit

264

sources is equal to the initial sea-surface height change and the displacement occurs

265

instantaneously at time t = 0 s. After the tsunami simulation, we calculate the pressure change p

266

at each OBPG location by subtracting the seafloor vertical movement from the sea-surface height267

change (Tsushima et al., 2012), ! = #! $! (& − (" ) (ρ0: seawater density, g0: gravitational

268

acceleration, η: sea-surface height change, and uz: seafloor vertical displacement). We here269

suppose ρ0 ~ 1.02 g/cm3 and g0 = 9.8 m/s2, so that 1 cm change of seawater column height (& −270(" ) equals to 1 hPa pressure change (i.e., ρ0g0 = 1 hPa/cm). We finally apply the same bandpass

271

filter to the simulated waveform as that applied to the observation.

manuscript submitted to Journal of Geophysical Research: Solid Earth

272

In the inversion analysis, we use the time-derivative waveforms of the bandpass-273

filtered pressure (!"/!$, the method of Kubota et al. (2018b)), because the time-derivative can274

reduce the artefacts due to the tsunami-irrelevant steps, which becomes the impulse and thus275

does not contain the offset change. The data time window used for the modeling is determined276

based on the visual inspection which includes the main part of the tsunami (listed in Table S1,

277

blue traces in Figure 3c). The goodness of the estimated source is evaluated using the variance278

reduction (VR):279(

VR = (1 −280∑! "#!"#$ $#!%&' %∑! #!"#$(+ × 100 (%)

(1)

281
282

where 2&'() and 2&*+, are the i-th observed and calculated data, respectively. We impose the283smoothing constraint for the inversion, and its weight is determined based on the trade-off

284

between the weight and the VR (Figure S1) to avoid both overfitting and oversmoothing.

285
286
287

4.2 Results
Figure 3 shows the results of the inversion. A subsidence with a horizontal extent of288~40 km × ~20 km, having a sharp peak near the GCMT centroid, was obtained (Figure 3a). The

289

direction of the northeast-southwest extents of the subsidence is consistent with the GCMT strike290

angle of 49°. The western edge of the subsidence region is consistent with the locations where291

the seafloor displacements of 1–2 m and fresh seafloor cracks were found by a seafloor

292

bathymetry survey just after the off-Fukushima earthquake conducted by Japan Agency for

293

Marine-Earth Science and Technology [JAMSTEC] (blue triangles in Figure 3a). The time294

derivatives of the S-net pressure waveforms were well reproduced (VR = 95.7%, Figure 3c).

295

Except for the steps just after the focal time at some near-source OBPGs, the observed pressure296

is also well explained (Figure 3b). The waveforms recorded at the other tsunami stations (Figure2971a) are also reproduced surprisingly well (Figures 3d–3g), even though they were not used for298

the inversion. This suggests that the use of the S-net data provides good spatial resolution of the299

tsunami source, and thus it is expected that we can obtain a reliable fault model. Note that the300

later arrivals in some stations (e.g., ~100 min at TM1 and TM2) are not well reproduced, which301

are caused by the coastal-reflections (Gusman et al., 2017). This is probably because the spatialmanuscript submitted to Journal of Geophysical Research: Solid Earth

302

resolution of the coastal shape from the topography data in our simulation is not sufficient to303

reproduce the reflected tsunami waves, and the high-resolution bathymetry data is necessary

304(Gusman et al., 2017; Kubota et al., 2018a).

305

306
307

Figure 3. Results of the tsunami source inversion. (a) Spatial distribution of the tsunami source.

308

Contour lines denote the subsided region with intervals of 20 cm. The white star is the JMA

309

epicenter, and blue triangles denote the location of the seafloor survey, where fresh surface

manuscript submitted to Journal of Geophysical Research: Solid Earth310

cracks were found. The yellow and gray circles show the S-net OBPGs used or not used,311

respectively, for inversion analysis. Comparisons of (b) the pressure waveforms and (c) the time-312

derivative waveforms. The gray and red traces denote the observed waveforms and simulated

313

waveforms from the tsunami source model. Traces marked by blue lines denote the time window314

used for the inversion analysis. (d–g) Waveform comparisons for the other networks, for (d)315

NOWPHAS Near-coastal GPS buoys, (e) NOWPHAS wave gauges, (f) OBPGs installed by ERI,

316

and (g) OBPGs installed by Tohoku University. At station 801, the waveforms during the data

317

missing are not drawn. See Figure 1 for station locations.

318
319

In the inversion, we used the time-derivative waveforms of the pressure to reduce the320

artefacts attributed to the tsunami-irrelevant pressure components (Kubota et al., 2018b). In order321

to see how well this method reduced the artefacts, we also conduct the additional inversion using322

the original pressure waveforms, instead of its time-derivative waveforms (Figure S2). The

323

weight of the smoothing is also determined based on the VR between the observed and simulated324

pressure waveforms (Figure S1). As a result, the distribution of the tsunami source is325

fundamentally similar to the original distribution (VR = 87 %), although a significant artificial

326

uplift of > 60 cm is estimated around S2N14 where the large step was recorded. Related to this327

artefact, the artificial pressure step irrelevant to the tsunami at S2N14 (see Section 3) was

328

reproduced, whereas that was not reproduced when using the time-derivative waveforms. In

329

order to avoid the artefact due to the tsunami-irrelevant step signals, the inversion using the time-330

derivative waveforms worked well.

331

We compare the tsunami source model estimated by the present study with the models

332

obtained using the tsunami data except for the S-net data (Figure S3, Table 1). The horizontal

333

location and spatial extent of the subsided region of our tsunami source model roughly334

correspond to those obtained by the previous studies. However, the amount of the maximum

335

subsidence was much larger than the previous models and the locations of the peak subsidence of336

tsunami source are slightly different from each other. Our tsunami source model had a maximum

337

subsidence of ~238 cm, whereas the two models obtained from the far-field tsunami data

338

(Adriano et al., 2018; Gusman et al. 2017) underestimated the subsidence (~180 cm and ~130339

cm, respectively, Table 1). The subsidence peak of our model was located ~10 km southeast and340

east of the models by Gusman et al. (2017) and Adriano et al. (2018), respectively (Figures S3a,

manuscript submitted to Journal of Geophysical Research: Solid Earth341

S3b). The peak location of the model of Nakata et al. (2019) was located ~10 km northwest of342

our models (Figure S3c), which was estimated by horizontally shifting the location of the slip

343

distribution model from the teleseismic data estimated by JMA (shown later, in Section 5.2) to fit

344

the coastal tide gauge waveforms. One reason for these differences may be the assumption of the345

fault geometry, but the more significant reason should be the station coverage and the source-346

station distance. The coastal tide gauges or the offshore stations used in these previous studies

347

were located far from the source region and the stations at the offshore side of the source region348

were not used in these studies, whereas the S-net has better station coverage and a smaller349

source-station distance. This could provide a better constraint on the horizontal location and peak

350

displacement amount to reproduce surprisingly well the tsunami waveforms not used for the351

inversion. Thanks to this improvement in the constraint, we believe that we can obtain a finite

352

fault model with a higher resolution, as shown in the next section.

353
354

5 Fault modeling355

5.1 Rectangular fault model with uniform slip

356

Here, we attempt to constrain the finite fault model of the off-Fukushima earthquake.

357

The horizontal location and the peak subsidence of our tsunami source distribution are slightly

358

different from the other previous models. Therefore, we first estimated the fault model based on359

a grid-search approach (Kubota et al., 2015; 2019). The procedure of the analysis is summarized360

in Text S2 but briefly explained here. We assume one planar rectangular fault with a uniform

361

slip, with a set of unknown parameters (the fault model candidate). The unknown parameters are362

the fault location (longitude, latitude, and depth) and its dimensions (length L and width W). The

363

strike, dip, and rake angles are fixed to the GCMT value (Table 1). The slip amount on fault D364

was adjusted to maximize the VR (Eq. (1)). The search range for the unknown parameters is

365

summarized in Table S2, determined based on the tsunami source model. Using the fault model

366

candidate, we calculated the seafloor displacement (Okada, 1992) and then calculated the367pressure changes at each OBPG. This calculation procedure is identical to that for the Green's368

function of the tsunami source modeling (see Text S1). The goodness of each fault model

369

candidate is evaluated by the VR values.

370
371

The horizontal location of the optimum fault is shown by dark red rectangle in Figure
4b. The detailed results are shown in Figure S4. We obtain the optimum fault model as L = 15

manuscript submitted to Journal of Geophysical Research: Solid Earth372

km, W = 10 km, and D = 467.7 cm (M0 = 2.1 × 1019 Nm, Mw 6.8, assuming a rigidity of μ = 30

373

GPa). The center of this model is located at a depth of 6 km, ~10 km east of the GCMT centroid

374(Table 1). The GCMT centroid depth was 12 km and the aftershocks are mainly located at depths

375of ~ 20 km (Figures 1b and 1c), whereas the estimated fault is located at the very shallow part of376

the crust (Figure 1c and Table 1). This disagreement has also been pointed out by Gusman et al.

377(2017), who suggested that the aftershocks determined from the inland network are378

systematically deeper than the actual depth. The horizontal extent of the tsunami source is379

relatively narrow and is located at the northeast, compared with the tsunami source model

380(Figure S4a). The reproductivity of the S-net pressure waveforms is reasonable (Figure S4),

381

although the VR is lower than that for the tsunami source modeling (VR = 59.3%). These382

mismatches are probably because of the simple assumption of the rectangular fault, which could383

not reproduce the southwest part of the tsunami source.

384

If we consider the empirical scaling relations from the magnitude, then the fault

385

dimension is expected to be ~700 km2 (e.g., Wells & Coppersmith, 1994). On the other hand, the386

estimated fault dimension of 150 km2 is much smaller. In order to assess the dimensions of the387

rectangular fault, fixing the seismic moment M0 and the fault center location to the optimum

388

model and varying the fault dimensions, we simulate tsunamis and compare the waveforms of389

representative S-net stations near the focal area (Figure S5). If we assume a larger fault with L >390

20 km, the arrival of the peak downheaval wave and its amplitude cannot be explained for the391

stations located northward (S2N01 and S2N02) or southward (S2N14 and S2N15) from the392source. In addition, the sharp peak of the downheaval waves observed at the stations located393

eastward (S2N09, S2N10, S2N11, S2N12, and S2N15) from the source are not well reproduced394

by the fault width for the case in which W > 15 km. These results suggest that the fault

395

dimensions should be L ≤ ~20 km and W ≤ ~15 km. Considering this range, the estimated fault

396

dimensions are obviously smaller than expected based on the scaling relation. These much

397

smaller fault dimensions are consistent with the size of the asperity, defined as the region of the398large slip on the fault (e.g., Somerville et al., 1999), expected from the empirical relation

399

deduced from the inland crustal earthquakes (Somerville et al., 1999; Miyakoshi et al., 2020).

400

This may suggest that this optimum rectangular fault corresponds to the asperity.

401
402

5.2 Slip distribution

manuscript submitted to Journal of Geophysical Research: Solid Earth403

We then conduct a finite fault inversion to estimate the slip distribution (finite fault

404

model). We assume a rectangular planar fault with dimensions of 45 km × 30 km, which passes

405

through the optimum fault obtained by the grid search, and then the planar fault is divided into

406

subfaults with size 3 km × 3 km. We then simulate the Green’s function, (the pressure change407

waveforms excited by each subfault), using a similar calculation procedure to that explained

408

above (see Texts S1 and S2). We estimated the slip amount of each subfault by the inversion

409

scheme identical to the tsunami source modeling, but we imposed a nonnegativity constraint

410(Lawson & Hanson, 1974). The smoothing constraint was also imposed (Figure S6). The other411

details are described in Text S3.

412

The slip distribution obtained by the inversion analysis and the tsunami source

413

distribution calculated from this slip distribution are shown in Figures 4a and 5a, respectively.

414

The tsunami source distribution (Figure 5a) is similar to that obtained by the tsunami source

415

inversion (Figure 3). The S-net and other tsunamis waveforms are explained (VR = 72.4%,

416

Figures 5b–5g). We obtain a maximum slip of Dmax = 637.2 cm, and the total seismic moment is417

M0 = 6.3 × 1019 Nm (Mw 7.1, μ = 30 GPa). The large slip is concentrated in the northeastern part418

of the fault plane, corresponding to the rectangular fault estimated by the grid-search analysis.

419

More specifically, subfaults with slip amounts with D > 0.5 × Dmax roughly correspond to the420

rectangular fault (subfaults marked by green lines in Figure 4a, 41% of the total M0, Mw 6.9). In

421

addition, a relatively small slip also extends to the southwestern part, which was not resolved in422

the grid-search analysis, probably because of the simple assumption of the uniform slip

423

rectangular fault. If we take subfaults with slip amounts larger than 0.2 × Dmax, both large

424

northeastern slip and relatively small southwestern slip are included (indicated by the thick black425

lines in Figure 4). Thus, we define these subfaults as the main rupture area. The main rupture

426

area had dimensions of ~30 km × ~ 20 km, and 81% of the total moment was concentrated in the427main rupture area. We calculate the centroid location xc = (xc, yc, zc) from the fault model, based428

on the slip-weighted average of subfault center locations over the main rupture area (pink star in

429

Figure 4a, Table 1), defined as:

430
431
432

3- =

∑ ! .! / !
∑ ! .!

,(2)manuscript submitted to Journal of Geophysical Research: Solid Earth433

where the subscript i denotes the subfault index, xi = (xi, yi, zi) is the center location and Di is the434

slip amount. We here used the subfaults within the main rupture area for the centroid calculation.

435

The horizontal location of the centroid is located ~5 km southeast from the GCMT centroid.

436

437
438

Figure 4. Result of the slip inversion. (a) Slip distribution (colored tiles). The pink and white439

stars indicate the slip-weighted averaged centroid and the JMA epicenter, respectively. Subfaults

440

with slip amounts larger than 0.2 × Dmax (the main rupture area) and larger than 0.5 × Dmax are

441

marked by thick black lines and green lines, respectively. (b) Shear stress change along the fault.

442

Negative (blue) and positive (red) values denote the shear stress decrease (or positive stress drop)443

and increase (negative stress drop), respectively. The dark red rectangle denotes the optimum444

rectangular fault obtained by the grid-search analysis.

445

manuscript submitted to Journal of Geophysical Research: Solid Earth

446
447

Figure 5. (a) Spatial distribution of the tsunami source calculated from the finite fault model. (b–448g) Comparisons of the observed and simulated waveforms. See Figures 3 and 4 for the other449

detailed description.

450
451

We then evaluated the resolution of the slip dimension from the S-net data based on the452

recovery test (Figure S7). In the recovery test, we simulated pressure waveforms assuming the453

rectangular fault with slip amount of 100 cm and with various fault dimensions (3 km, 6 km, 9manuscript submitted to Journal of Geophysical Research: Solid Earth454

km, 12 km, and 15 km). Then, these simulated waveforms are regarded as the observed data and455

inverted to estimate the slip distribution. The other conditions are identical to the original456

inversion. As a result, when assuming the smaller fault dimension (3 or 6 km), the estimated457

maximum slip was much smaller than the input and the recovered slip spreads out around the458

surroundings (Figure S7a–S7b). On the other hand, when assuming the faults with dimension of459≥ ~9 km, the slip smearing is small and the recovered maximum slip is consistent with the input

460

amount (Figure S7c–S7f). This suggests that the resolution of the slip distribution is ~9 km.

461

We compare our fault model with the other models estimated using the teleseismic data462

and onshore geodetic data, by JMA and Geospatial Information Authority, Japan (GSI),

463

respectively (Figure S8). The detailed fault parameters of these models are summarized in Table

464

1. The center of these fault models are ~10 km, almost consistent but slightly deeper compared465

with our model (~8 km). The location of the maximum slip in the JMA teleseismic fault model is466

almost identical to our model, although the maximum slip amount of 4.0 m was much smaller

467

than ours (6.4 m) (Figure S8a). The fault length of the geodetic fault model by GSI was ~45 km,

468

which was larger than the length of the main rupture area of our model (~30 km). The maximum

469

slip of GSI model, 0.78 m, was much smaller than the average slip amount within the main

470

rupture area of our model (2.8 m) (Figure S8b). This comparison possibly suggests the S-net data471

improved the spatial resolution of the slip distribution.

472

We then simulated the pressure waveforms for the representative S-net stations using the473

previously-proposed fault models (Adriano et al., 2018; Gusman et al., 2017; Nakata et al., 2019,

474

Figure S9). We found the waveforms at the northern stations (~38°N) far from the focal area are475

reasonably reproduced by these model (Figures S9a) although the arrival times at the waveforms

476

near the focal area (~36.5–37°N) and at the south (~36°N) were not reproduced (Figures S9b–477

S9c). The impulsive short-wavelength tsunami features at stations to the south direction (S2N15,

478

S2N26, and S1N01) were not also explained well.

479

To investigate how the S-net data improved the resolution of the finite fault model, we

480then conducted the additional finite fault inversion using only the stations far from the focal area

481

(Figures S10 and S11). The waveforms used for this additional inversion analysis is marked by482

blue traces in Figure S11c. The other settings is identical to the original analysis. As a result, the483

maximum slip of 3.8 m is almost similar to that obtained in the past studies (Figure S10a, S10c–484

S10e), but the large slip peak at the northeast part of the fault, around the rectangular fault

manuscript submitted to Journal of Geophysical Research: Solid Earth485

estimated by the grid-search analysis (Section 5.1), was not resolved (Figure S10a). The

486

waveforms for the far-field stations used for the inversion were reasonably explained, although487

the short-period impulsive tsunamis at the near-field S-net stations were not reproduced (Figure488

S11). In addition, we also conducted the recovery test to evaluate the resolution (Figure S12),

489

with the identical procedure to that in the original analysis (Figure S7). We found that the490

recovered slip image was unsharp and the resolution was much lower than the original (> ~15

491km). This result indicates the near-source S-net OBPGs are important to resolve the slip

492

distribution and large slip component at the asperity with high resolution.

493

In our fault model, the downdip depth of the main rupture area was estimated as12.9–

494

14.6 km. We evaluated the resolution of the downdip limit of the fault by additional tsunami

495

simulations (Figure S13). We assume a simple rectangular fault with the length of 30 km and slip

496

amount of 3 m, based on the main rupture area of the finite fault model (Figure S13a). We then497

vary the fault width so that the downdip limit of the fault is varied, and simulate tsunamis. As a498

result, the waveforms at the southern stations are different depending on the fault width,499

particularly for the stations located at the southeast from the epicenter (e.g., S2N18, S2N19,

500

Figures S13c–S13d). At these stations, the simulated first down waves arrive earlier and its

501

duration is longer, when the larger fault width are assumed (corresponding fault bottom depth is502≥ ~16 km), whereas the down waves are reproduced when the fault bottom end is assumed as503

13–15 km. This suggests that the down dip end of the fault should be ≤ 15 km, consistent with504

the depth obtained in the finite fault modeling. However, the simulated waveforms are almost

505

similar regardless of the fault width at the stations at the north (e.g., S3N26, S2N01, Figure

506S13b). Taking the point into account that the previous fault models (Figures S10c–S10e) were507

derived without using the stations located at the offshore side of the focal area, the use of the508

high-coverage S-net data, particularly located at the southeast of the focal area, potentially509

contributed to the constraint of the downdip depth of the fault of the off-Fukushima earthquake.

510
511
512

5.3 Stress drop
In Figure 4b, we calculate the distribution of the shear stress change along the fault (i.e.,

513

stress drop) by computing the shear stress change at the center of each subfault using the514

equation of Okada (1992). Here the shear stress change is calculated on fault plane along the slip

515

vector (strike = 49°, dip = 36°, rake = −89°). The rectangular fault estimated by the grid-searchmanuscript submitted to Journal of Geophysical Research: Solid Earth516

analysis agrees well with the region where the shear stress is largely released (dark red rectangle

517

in Figure 4b), indicating that the rectangular fault model corresponds to the asperity, as discussed

518

above. We then calculate the energy-based stress drop, or the slip-weighted average stress drop,519

ΔσE (Noda et al., 2013) as:

520
521

Δ40 =

∑! .! Δ1!
∑ ! .!

,

(3)

522
523

where Di is the slip amount at the i-th subfault, and Δσi is the stress drop at the i-th fault. Using524

the subfaults within the main rupture area (D > 0.2 × Dmax), we obtain ΔσE = 4.2 MPa. As it is a525

difficult issue to choose the appropriate area for the ΔσE calculation (Brown et al., 2015), we also526

calculate the stress drop using all subfaults and the subfaults with D > 0.5 × Dmax (the region

527

marked by green lines in Figure 4a), and we obtain ΔσE = 3.3 MPa and ΔσE = 6.8 MPa,

528

respectively. This value seems not so small as expected for the interplate earthquakes (~100 MPa,529

e.g., Kanamori & Anderson 1975), but rather is consistent with the intraplate earthquakes, which

530

generally have larger stress drop values (e.g., Miyakoshi et al., 2020; Somerville et al., 1999).

531

In Figure S10b, we calculate the spatial distribution of the shear stress change from532

based on the fault model estimated by the additional inversion which uses only the far-field

533

stations (Figure S10a). Using the subfaults within the area where the slip was larger than 20 % of534

the maximum slip, we obtain the average stress drop of ΔσE = 2.3 MPa. However, because the S-

535

net OBPGs near the source were not explained by this additional fault model (Figure S11), the536

actual average stress drop value should be larger than this value.

537
538
539

6. Discussion: implication for the intraplate stress regime
After the 2011 Tohoku earthquake, it has been reported that the normal-faulting

540

seismicity significantly increased in the upper plate, which is thought to be related to the stress541

perturbation by the Tohoku earthquake (Figures 1d–1f, Asano et al., 2011; Hasegawa et al.,542

2012; Yoshida et al., 2012). This change in seismicity is interpreted as the result whereby the543

intraplate stress regime switched after the Tohoku earthquake from the horizontal compression to544

the horizontal extension (e.g., Hasegawa et al. 2012). As discussed previously, the use of the S-

545

net tsunami data improved the constraint on the tsunami source and the fault model of the off-

manuscript submitted to Journal of Geophysical Research: Solid Earth546

Fukushima earthquake, which made it possible to obtain the detailed distribution of the shear

547

stress reduction and the static stress drop. Using these results, we attempt to discuss the548quantitative relationship between the crustal stress released during the off-Fukushima earthquake

549

and the stress increase due to the 2011 Tohoku earthquake. This kind of discussion is typically550

difficult to conduct because it is rare that both the high-resolution fault model of the M~7551

offshore earthquake and the significant stress perturbation due to the megathrust earthquake are552

available.

553

If the stress regime switched from the horizontal compression to the horizontal

554

extension by the Tohoku earthquake around the off-Fukushima earthquake, the deviatoric stress,555

or the initial shear stress on the fault of the off-Fukushima earthquake, should be smaller than (or556

at least equivalent to) the static shear stress increase due to the Tohoku earthquake (Figure 6a).

557

In other words, the stress drop of the off-Fukushima earthquake should be smaller than the shear

558

stress increase due to the Tohoku earthquake. In Figure 7a, we calculate the shear stress change559

due to the Tohoku earthquake, using the fault model of Iinuma et al. (2012), along the fault

560

geometry of the off-Fukushima earthquake. The shear stress change related to the Tohoku

561

earthquake around the focal area of the off-Fukushima earthquake is only ~2 MPa, which is562

smaller than the stress drop of the off-Fukushima earthquake. We also calculated the shear stress563

change related to the Tohoku earthquake based on the fault model of Yamazaki et al. (2018), as564~1.1 MPa. In addition, at the some subfaults in which the large slip was estimated (e.g., the565

region marked by green lines in Figure 4a, D > 0.5 × Dmax), the stress drop values are higher (>566~10 MPa, the slip-weighted average is ΔσE = 6.8 MPa) and thus the stress discrepancy was much

567

larger. The larger stress drop of the off-Fukushima earthquake than the stress increase after the568

Tohoku earthquake is inconsistent with the presumption that the intraplate stress regime switched569

by the stress change of the Tohoku earthquake. There may be other causes for the normal-

570

faulting stress regime around the focal area, particularly at the asperity where the larger slip was

571

estimated, of the off-Fukushima earthquake.

572

manuscript submitted to Journal of Geophysical Research: Solid Earth

573
574

Figure 6. Schematic illustration of the temporal change of the stress regime around the off-575

Fukushima earthquake. The stress regimes (a) assuming the switching of the stress regime after

576

the Tohoku earthquake and (b) without assuming the stress switching.

577

578
579

Figure 7. Horizontal distribution of the shear stress change on the 2016 off-Fukushima

580

earthquake fault at a depth of 10 km. Shear stress changes along the fault geometry of the off-581

Fukushima earthquake due to (a) the Tohoku earthquake (Iinuma et al. 2012), (b) postseismic582

slip during April 2011 and December 2011 (Iinuma et al., 2016), and (c) Postseismic slip during

583

September 2012 and September 2016 (Tomita et al. 2020). The main rupture area of the 2016

manuscript submitted to Journal of Geophysical Research: Solid Earth584

off-Fukushima earthquake inferred from the inversion analysis is also indicated by black lines.

585

Note that the color scales are different in each subfigure.

586
587

One possible cause is the postseismic slip of the Tohoku earthquake (Iinuma et al.,588

2016; Tomita et al., 2020). After the Tohoku earthquake, the postseismic seafloor deformation589

was detected by the seafloor geodetic observation (Tomita et al., 2015; 2017), which was caused590

by the postseismic slip along the fault and the viscoelastic deformation (Iinuma et al., 2016; Sun591

et al., 2014; Tomita et al., 2020). Among the postseismic deformation, the afterslip along the592

plate interface is dominant in the south of the rupture area of the Tohoku earthquake, including593

the region off Fukushima, whereas the viscoelastic deformation dominates the northern part of594

the Tohoku earthquake rupture area (Iinuma et al., 2016; Tomita et al., 2020). We calculate the595

shear stress change on the 2016 fault geometry using the postseismic slip models to evaluate the596

contribution by the postseismic slip around the focal area. We calculate the stress change due to597

the postseismic slip models from 23 April 2011 to 10 December 2011 (Iinuma et al., 2016,

598

Figure 7b) and during 2012 and 2016 (Tomita et al., 2020, Figure 7c), but the contributions by599

these postseismic slip models were minor (on the order of 10−1 MPa). We therefore concluded

600

that the shear stress increase due to the postseismic slip could not resolve the apparent

601

contradiction between the stress drop of the off-Fukushima earthquake and the shear stress602

increase after the Tohoku earthquake. This contradiction may arise from the assumption of the603

switching of the stress regime, which was a reverse-faulting and a normal-faulting regime before

604

and after the Tohoku earthquake, respectively.

605

It will be reasonable to interpret this apparent contradiction that the horizontal

606

extensional stress regime was already predominant around the 2016 off-Fukushima earthquake

607

even before the Tohoku earthquake and the stress increase by the Tohoku earthquake further

608

enhanced the extensional stress (Figure 6b), in contrary to the past studies which report the609

horizontal compressive stress attributed to the plate coupling force was widely dominant in Japan610

before the Tohoku earthquake (e.g., Terakawa & Matsu’ura, 2010; Wang & Suyehiro, 1999). In

611

other words, the extensional stress accumulated even before the 2011 Tohoku earthquake was the612

dominant cause for the 2016 off-Fukushima earthquake. If we assume the stress state within the613

plate switched from the horizontal compression to the extension due to the stress change of a few

614

MPa by Tohoku earthquake (Figure 6a), the 2016 earthquake must have occurred under the verymanuscript submitted to Journal of Geophysical Research: Solid Earth615

low extensional stress level (less than a few MPa), but this is very unlikely to occur. Although616

there will be some uncertainties in the stress drop estimation of the 2016 earthquake and the617

stress increase due to the 2011 Tohoku earthquake, these uncertainties do not matter to this618

discussion and it is very reasonable to interpret that the stress state before the 2011 Tohoku

619

earthquake was the horizontal extensional regime (Figure 6b). There are some recent reports that620

some normal-faulting microearthquakes occurred even before the Tohoku earthquake in the621

inland region of Fukushima prefecture (Imanishi et al., 2012; Yoshida et al., 2015a; 2015b),622

which supports our hypothesis of the normal-faulting stress regime being predominant in this623

location even before the Tohoku earthquake.

624

One possible reason for this normal-faulting stress regime is the effect of bending of625

the overriding plate, in which the horizontal extensional and compressional stresses develop at626

the shallower and the deeper portion of the plate, respectively (e.g., Fukahata & Matsu’ura,

627

2016; Hashimoto & Matsu’ura 2006; Turcotte & Schubert, 2002). Yoshida et al. (2015a) showed628that the normal-faulting stress regime is dominant at depths shallower than ~15 km in this region,629

while the reverse-faulting stress regime is dominant at depths greater than ~15 km, which is630

consistent with the hypothesis. We can also consider the topographic effects (Sasajima et al.,631

2019; Wang et al., 2019) for the formation of the horizontal extensional stress. There may be

632

another possible interpretation for this contradiction, that the stress regime switched to the633

reverse-faulting regime again by the off-Fukushima earthquake; However this is improbable634

because normal-faulting seismicity can be found nearby, even after one year from the earthquake

635(Figure 1f).

636

Some major normal-faulting earthquakes were reported around the focal area of the637

off-Fukushima earthquake in 1938 (Abe, 1977; Murotani, 2018). Furthermore, according to the638

geologic cross-section around the off-Iwaki gas field, located near the 2016 off-Fukushima

639

earthquake, the northeast-southwest-trending reverse faults were developed at a depth shallower

640

than 6 km, which are considered to have formed during Oligocene and Miocene (Iwata et al.,

641

2002). Along this fault trace, it was also reported that the normal-faulting-type surface offsets

642with vertical offset of 5–10 m were found, and it was suggested that the direction of the tectonic

643

stress flipped to the normal-faulting regime during Quaternary and normal-faulting earthquakes

644

similar to the 2016 off-Fukushima earthquake repeatedly occurred along this fault (S. Toda,645

https://irides.tohoku.ac.jp/media/files/earthquake/eq/2016_fukushima_eq/20161122_fukushima_manuscript submitted to Journal of Geophysical Research: Solid Earth646

eq_activefault_toda.pdf, in Japanese). These reports may support our hypothesis that the crustal

647

stress regime was under the normal-faulting regime even before the Tohoku earthquake.

648

Note that the downdip limit of the main rupture area of our fault model the off-649

Fukushima earthquake is estimated as < ~15 km (Figure S9), which is approximately consistent

650

with the downdip limit depth of the normal-faulting regime in the inland Fukushima region651

estimated by Yoshida et al. (2015a). This suggests that the horizontal extensional stress regime652

before the Tohoku earthquake around the focal area of the off-Fukushima earthquake is653

predominant at depths shallower than 15 km and the stress neutral zone related to bending of the654

overriding plate lies at a depth of ~15 km. However, we note that the normal-faulting seismicity

655

extensively increased in the overriding plate after the Tohoku-Oki earthquake even at a depth

656

deeper than 15 km (e.g., Asano et al., 2011; Hasegawa et al., 2012). This might suggest that the657

stress-neutral depth deepened around this region after the Tohoku earthquake.

658

As a summary of this discussion, the temporal change of the intraplate crustal stress659

around the off-Fukushima earthquake can be interpreted as follows. The horizontal extensional

660

stress was predominant before the Tohoku earthquake within the shallowest part of the661

continental plate, but may not exceed the crustal strength. After the Tohoku earthquake, its stress662

perturbation enhanced the extensional stress, provoking the normal-faulting seismicity.

663

Before the 2011 Tohoku earthquake, no major seismicity was detected around the focal664

area of the off-Fukushima earthquake (e.g., Asano et al., 2011; Hasegawa et al., 2012) and the665

onshore seismic network could not detect micro-seismicity around this offshore region. On the666

other hand, the use of the S-net OBPGs could well constrain the fault modeling of the 2016 off-667

Fukushima earthquake, which provides an important implication for the crustal stress regime668

prior to the Tohoku earthquake, even though the S-net was not installed at that time. Such669

information about the stress regime is important to understand the spatio-temporal change of the670

intraplate stress state and the generation mechanisms of the intraplate earthquake, especially after

671

a megathrust earthquake. Our analysis demonstrated that the analysis of the offshore S-net data

672

provided implications for the crustal stress regime at the offshore region, which was difficult to673

discuss before the S-net was available. Although the S-net OBPG data contains the tsunami-

674

irrelevant pressure change signals, careful analysis of this data significantly improves the675

constraint of the fault model and will deepen our understanding of the earthquake generation.

676

manuscript submitted to Journal of Geophysical Research: Solid Earth

677

7 Conclusions

678

We examined the S-net tsunami data associated with the off-Fukushima earthquake on

679

21 November 2016 (Mw 7.1). We first processed the S-net OBPG data and found some pressure

680

signals irrelevant to tsunami were observed: (1) an extremely large drift component and (2) an

681

abrupt pressure step around the origin time. We discussed the cause of these tsunami-irrelevant

682

signals and concluded that these signals were not due to the pressure sensors themselves but

683

probably due to the observation system, although further investigations are necessary. We then684

analyzed the S-net data in order to estimate the tsunami source model and the fault model.

685

Careful analysis of the S-net OBPG data provided the tsunami source distribution, which had a686

large subsidence with strike angle consistent with the GCMT solution. Our fault model suggested

687

that the energy-based stress drop of the off-Fukushima earthquake is ΔσE ~4.2 MPa. The

688

quantitative comparison between the stress drop and the static stress changes caused by the 2011

689

Tohoku earthquake and its postseismic slip suggested that the additional source of the horizontal

690

extensional stress is necessary to explain the stress drop. We interpreted the stress regime around691

the off-Fukushima earthquake to be the horizontal extensional even before the Tohoku

692

earthquake, related to the bending of the overriding plate. The S-net pressure data is very useful693

to constrain the tsunami source model and the finite fault model, even if the model is perturbed694

by the tsunami-irrelevant signals, which provided an important implication for the tectonic stress695

regime within the overriding plate.

696
697
698

Data Availability Statement
The S-net OBPG data (National Research Institute for Earth Science and Disaster

699

Resilience [NIED], 2019, https://doi.org/10.17598/NIED.0007) are available with data request

700

and permission, through the website of NIED (https://www.seafloor.bosai.go.jp, in Japanese).

701

The data policy of the S-net data is shown in702

https://www.mowlas.bosai.go.jp/policy/?LANG=en. The NOWPHAS tsunami data is provided

703

upon request to the Port and Airport Research Institute (PARI), in which the data redistribution is704

prohibited (the contact address is shown in https://nowphas.mlit.go.jp/pastdata/, in Japanese,

705

accessed on 21 July, 2021). The data of the OBPGs installed by ERI of the University of Tokyo,

706

as used in Gusman et al. (2017), are attached to this article as Supplementary Dataset S1, with707

permission of Masanao Shinohara, ERI. The OBPG data of Tohoku University are available inmanuscript submitted to Journal of Geophysical Research: Solid Earth708

Hino et al. (2021, https://doi.org/ 10.5281/zenodo.4961355). Station locations of the S-net OBPG709

are available at https://www.seafloor.bosai.go.jp/st_info/. The location of the OBPGs installed by

710

the ERI is available in Gusman et al. (2017), as well as in Dataset S1. The locations of the711

NOWPHAS GPS buoys and wave gauges are available at https://nowphas.mlit.go.jp/pastdata/ (in712Japanese). The locations of the OBPGs installed by Tohoku University are listed in Hino et al.

713(2021), as well as in Table S3.

714

We purchased the JTOPO30v2 bathymetry data from the Marine Information

715

Research Center (http://www.mirc.jha.jp/en/) of the Japan Hydrographic Association. The plate716

boundary model in Figure 1 (Nakajima & Hasegawa, 2006) is available from the website of717

Fuyuki Hirose (https://www.mri-jma.go.jp/Dep/sei/fhirose/plate/PlateData.html, in Japanese,718

accessed on 21 July, 2021). The rotation data of the S-net sensor (Takagi et al., 2019) was

719

provided by contacting Ryota Takagi, the lead author of Takagi et al. (2019). The slip models of720

the mainshock and postseismic slip of Iinuma et al. (2012; 2016) and Tomita et al. (2020) were721

provided by the corresponding authors of each article, Takeshi Iinuma and Fumiaki Tomita. The

722

slip distribution models of Gusman et al. (2017), Adriano et al. (2018), and Nakata et al. (2019)723

are available in each paper. The finite fault model using the teleseismic data by JMA is shown in724

https://www.data.jma.go.jp/svd/eqev/data/sourceprocess/event/2016112205594689far.pdf (in725Japanese, accessed on 21 July, 2021) and its digital data is available from726

https://www.data.jma.go.jp/svd/eqev/data/sourceprocess/data/2016112205594689far.zip

727(accessed on 21 July, 2021). The fault model using the onshore geodetic data by GSI is available

728

from https://cais.gsi.go.jp/YOCHIREN/activity/214/214.e.html (accessed on 21 July, 2021) and729

https://cais.gsi.go.jp/YOCHIREN/activity/214/image214/008.pdf (in Japanese, accessed on 21

730

July 2021). The location of the seafloor bathymetry survey conducted by the Japan Agency for

731

Marine-Earth Science and Technology [JAMSTEC] (blue triangles in Figure 4a) was taken from732

http://www.jamstec.go.jp/ceat/j/topics/20161208.html,

733

http://www.jamstec.go.jp/j/about/press_release/20170301/ (in Japanese, accessed on 21 July,7342021).

735

The results of the tsunami source modeling (Figure 3) and the slip distribution (Figure7364) are available in the tgz compressed file, Dataset S2. The detailed description of the dataset is737

available in README file of the tgz file.

738

manuscript submitted to Journal of Geophysical Research: Solid Earth739

Acknowledgments

740

The present study was supported by JSPS KAKENHI Grant Number JP19K14818 and741

JP19H02409 from the Japan Society for the Promotion of Science. Discussions with Ryota Hino,

742

Naoki Uchida, and Tatsuhiko Saito, about interpreting the S-net pressure signals and the results743

of the fault modeling, were fruitful. We would also like to thank Ryota Takagi, Takeshi Iinuma,

744

and Fumiaki Tomita for providing their useful data. We thank the Editor Rachel Abercrombie

745

and the Associate Editor, and two anonymous reviewers for their comments.

746

manuscript submitted to Journal of Geophysical Research: Solid Earth747

References

748

Abe, K. (1977). Tectonic implications of the large shioya-oki earthquakes of 1938. Tectonophys,749

41, 269–289. https://doi.org/10.1016/0040-1951(77)90136-6

750

Adriano, B., Fujii, Y., & Koshimura, S. (2018). Tsunami source and inundation features around751

Sendai Coast, Japan, due to the November 22, 2016 Mw 6.9 Fukushima earthquake.

752

Geoscience Letters, 5, 2. https://doi.org/10.1186/s40562-017-0100-9753

Aoi, S., Suzuki, W., Chikasada, N. Y., Miyoshi, T., Arikawa, T., & Seki, K. (2019).

754

Development and utilization of real-time tsunami inundation forecast system using s-

755

net data. Journal of Disaster Research, 14(2), 212–224.

756

https://doi.org/10.20965/jdr.2019.p0212

757

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

758

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

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

760

Asano, Y., Saito, T., Ito, Y., Shiomi, K., Hirose, H., Matsumoto, T., … Sekiguchi, S. (2011).

761

Spatial distribution and focal mechanisms of aftershocks of the 2011 off the Pacific762

coast of Tohoku Earthquake. Earth, Planets and Space, 63, 669–673.

763

https://doi.org/10.5047/eps.2011.06.016

764

Baba, T., Hirata, K., Hori, T., & Sakaguchi, H. (2006). Offshore geodetic data conducive to the765

estimation of the afterslip distribution following the 2003 Tokachi-oki earthquake.

766

Earth and Planetary Science Letters, 241, 281–292.

767

https://doi.org/10.1016/j.epsl.2005.10.019

768

Brown, L., Wang, K., & Sun, T. (2015). Static stress drop in the Mw 9 Tohoku-oki earthquake:

769

Heterogeneous distribution and low average value. Geophysical Research Letters,

770

42, 10,595–10,600. https://doi.org/10.1002/2015GL066361

771

Chadwick, W. W. J., Nooner, S. L., Zumberge, M. A., Embley, R. W., & Fox, C. G. (2006).

772

Vertical deformation monitoring at Axial Seamount since its 1998 eruption using773

deep-sea pressure sensors. Journal of Volcanology and Geothermal Research, 150,

774

313–327. https://doi.org/10.1016/j.jvolgeores.2005.07.006

775

Dhakal, Y. P., Kunugi, T., Suzuki, W., Kimura, T., Morikawa, N., & Aoi, S (2021). Strong776

motions on land and ocean bottom: Comparison of horizontal PGA, PGV, and 5%

777

damped acceleration response spectra in northeast Japan and the Japan Trench area.

manuscript submitted to Journal of Geophysical Research: Solid Earth778

Bulletin of the Seismological Society of America, 111.

779

https://doi.org/10.1785/0120200368

780

Fukahata, Y., & Matsu’ura, M. (2016). Deformation of island-arc lithosphere due to steady plate781

subduction. Geophysical Journal International, 204, 825–840.

782

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

783

Fukuyama, E., Ishida, M., Dreger, D. S., & Kawai, H. (1998). Automated seismic moment tensor

784

determination by using on-line broadband seismic waveforms. Zisin 2, 51, 149–156

785(in Japanese with English abstract). https://doi.org/10.4294/zisin1948.51.1_149

786

Gusman, A. R., Satake, K., Shinohara, M., Sakai, S., & Tanioka, Y. (2017). Fault slip

787

distribution of the 2016 Fukushima earthquake estimated from tsunami waveforms.

788

Pure and Applied Geophysics, 174, 2925–2943. https://doi.org/10.1007/s00024-017-

789

1590-2790

Hardebeck, J. L., & Okada, T. (2018). Temporal stress changes caused by earthquakes: a review.

791

Journal of Geophysical Research: Solid Earth, 123, 1350–1365.

792

https://doi.org/10.1002/2017JB014617

793

Hasegawa, A., Yoshida, K., Asano, Y., Okada, T., Iinuma, T., & Ito, Y. (2012). Change in stress794

field after the 2011 great Tohoku-Oki earthquake. Earth and Planetary Science

795

Letters, 355–356, 231–243. https://doi.org/10.1016/j.epsl.2012.08.042

796

Hashimoto, C., & Matsu’ura, M. (2006). 3-D simulation of tectonic loading at convergent plate797

boundary zones: Internal stress fields in northeast Japan. Pure and Applied

798

Geophysics, 163(9), 1803–1817. https://doi.org/10.1007/s00024-006-0098-y

799

Hino, R. (2015). An overview of the Mw 9, 11 March 2011, Tohoku earthquake. Summary of the800

Bulletin of the International Seismological Centre, 48, 100–132.

801

https://doi.org/10.5281/zenodo.998789

802

Hino, R., Suzuki, S., Sato, M., & Kubota, T. (2021) Ocean bottom pressure data of the 2016 Off-803

Fukushima earthquake obtained by Tohoku University [Data set]. Zenodo.

804

http://doi.org/10.5281/zenodo.4961355

805

Hino, R., Inazu, D., Ohta, Y., Ito, Y., Suzuki, S., Iinuma, T., … Kaneda, Y. (2014). Was the8062011 Tohoku-Oki earthquake preceded by aseismic preslip? Examination of seafloor

807

vertical deformation data near the epicenter. Marine Geophysical Research, 35, 181–808

190. https://doi.org/10.1007/s11001-013-9208-2

manuscript submitted to Journal of Geophysical Research: Solid Earth809

Hua, Y., Zhao, D., Toyokuni, G., & Xu, Y. (2020). Tomography of the source zone of the great

810

2011 Tohoku earthquake. Nature Communications, 11, 1163.

811

https://doi.org/10.1038/s41467-020-14745-8

812

Iinuma, T., Hino, R., Kido, M., Inazu, D., Osada, Y., Ito, Y., … Miura, S. (2012). Coseismic slip

813

distribution of the 2011 off the Pacific Coast of Tohoku Earthquake (M9.0) refined814

by means of seafloor geodetic data. Journal of Geophysical Research, 117, B07409.

815

https://doi.org/10.1029/2012JB009186

816

Iinuma, T., Hino, R., Uchida, N., Nakamura, W., Kido, M., Osada, Y., & Miura, S. (2016).

817

Seafloor observations indicate spatial separation of coseismic and postseismic slips

818

in the 2011 Tohoku earthquake. Nature Communications, 7, 13506.

819

https://doi.org/10.1038/ncomms13506

820

Imanishi, K., Ando, R., & Kuwahara, Y. (2012). Unusual shallow normal-faulting earthquake

821

sequence in compressional northeast Japan activated after the 2011 off the Pacific822

coast of Tohoku earthquake. Geophysical Research Letters, 39(9), L09306.

823

https://doi.org/10.1029/2012GL051491

824

Inazu, D., & Hino, R. (2011). Temperature correction and usefulness of ocean bottom pressure

825

data from cabled seafloor observatories around Japan for analyses of tsunamis, ocean

826

tides, and low-frequency geophysical phenomena. Earth, Planets and Space, 63(11),8271133–1149. https://doi.org/10.5047/eps.2011.07.014

828

Inoue, M., Tanioka, Y., & Yamanaka, Y. (2019). Method for near-real time estimation of829

Tsunami sources using ocean bottom pressure sensor network (S-net). Geosciences,830

9, 310. https://doi.org/10.3390/geosciences9070310

831

Iwata, T., Hirai, A., Inaba, T., & Hirano, M. (2002). Petroleum system in the Offshore Joban

832

Basin, northeast Japan. Journal of the Japanese Association for Petroleum

833

Technology, 67(1), 62–71 (in Japanese with English abstract).

834

https://doi.org/10.3720/japt.67.1_62

835
836

Kanamori, H., & Anderson, D. L. (1975). Theoretical basis of some empirical relations inseismology. Bulletin of the Seismological Society of America, 65(5), 1073–1095.

837

Kanazawa, T., & Hasegawa, A. (1997). Ocean-bottom observatory for earthquakes and tsunami

838

off Sanriku, north-eastern Japan using submarine cable. Proceedings of International

839

Workshop on Scientific Use of Submarine Cables, 208–209.

manuscript submitted to Journal of Geophysical Research: Solid Earth

840

Kanazawa, T., Uehira, K., Mochizuki, M., Shinbo, T., Fujimoto, H., Noguchi, S., et al. (2016). S‐

841

net project, cabled observation network for earthquake and tsunamis [Conference842

presentation]. Suboptic2016, Dubai. https://suboptic.org/wp‐content/uploads/fromke-843

vin/program/WE2B.3%20S‐

844

net%20Project,%20Cabled%20Observation%20Network%20for%20Earthquakes%2

845

0and%20Tsunamis.pdf

846

Kawaguchi, K., Sakuraba, S., Nagai, T., Nakai, K., Nukada, K., & Murase, H. (2017).

847

Characteristics of observed low frequency wave profiles of the 2016 Fukushima-Off

848

earthquake at offshore, shallow sea and coastal stations. Journal of Japan Society of849

Civil Engineers, Ser. B3 (Ocean Engineering), 73(2), I_743-I_748.

850

https://doi.org/10.2208/jscejoe.73.I_743

851

Kodaira, S., Fujiwara, T., Fujie, G., Nakamura, Y., & Kanamatsu, T. (2020). Large coseismic852

slip to the trench during the 2011 Tohoku-Oki earthquake. Annual Review of Earth853

and Planetary Sciences, 48, 321–343. https://doi.org/10.1146/annurev-earth-071719-

854

055216855

Kodaira, S., Iinuma, T., & Imai, K. (2021). Investigating a tsunamigenic megathrust earthquake856

in the Japan Trench. Science, 371, eabe1169. https://doi.org/10.1126/science.abe1169

857

Kubota, T., Hino, R., Inazu, D., Ito, Y., & Iinuma, T. (2015). Complicated rupture process of the858

Mw 7.0 intraslab strike-slip earthquake in the Tohoku region on 10 July 2011859

revealed by near-field pressure records. Geophysical Research Letters, 42, 9733–

860

9739. https://doi.org/10.1002/2015GL066101

861

Kubota, T., Saito, T., Ito, Y., Kaneko, Y., Wallace, L. M., Suzuki, S., Hino, R., & Henrys, S.

862(2018a). Using tsunami waves reflected at the coast to improve offshore earthquake863

source parameters: application to the 2016 Mw 7.1 Te Araroa earthquake, New

864

Zealand. Journal of Geophysical Research: Solid Earth, 123, 8767–8779.

865

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

866

Kubota, T., Suzuki, W., Nakamura, T., Chikasada, N. Y., Aoi, S., Takahashi, N., & Hino, R.

867(2018b). Tsunami source inversion using time-derivative waveform of offshore868

pressure records to reduce effects of non-tsunami components. Geophysical Journal

869

International, 215, 1200–1214. https://doi.org/10.1093/gji/ggy345

870

manuscript submitted to Journal of Geophysical Research: Solid Earth

871

Kubota, T., Hino, R., Inazu, D., & Suzuki, S. (2019). Fault model of the 2012 doublet

872

earthquake, near the up-dip end of the 2011 Tohoku-Oki earthquake , based on a873

near-field tsunami : implications for intraplate stress state. Progress in Earth and874

Planetary Science, 6, 67. https://doi.org/10.1186/s40645-019-0313-y

875

Kubota, T., Saito, T., & Suzuki, W. (2020a). Millimeter‐scale tsunami detected by a wide and876

dense observation array in the deep ocean: Fault modeling of an Mw 6.0 interplate877

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

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

879

Kubota, T., Saito, T., Chikasada, N. Y., & Suzuki, W. (2020b). Ultrabroadband seismic and880

tsunami wave observation of high‐sampling ocean‐bottom pressure gauge covering881

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

882

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

883

Kubota, T., Saito, T., Chikasada, N. Y., & Sandanbata, O. (2021). Meteotsunami observed by the884

deep-ocean seafloor pressure gauge network off northeastern Japan. [Preprint]. Earth885

and Space Science Open Archive. https://doi.org/10.1002/essoar.10506159.1

886

Lay, T. (2018). A review of the rupture characteristics of the 2011 Tohoku-oki Mw 9.1

887
888

earthquake. Tectonophysics, 733, 4–36. https://doi.org/10.1016/j.tecto.2017.09.022
Matsubara, M., Sato, H., Uehira, K., Mochizuki, M., Kanazawa, T., Takahashi, N., Suzuki, K., &889

Kamiya, S. (2019). Seismic velocity structure in and around the Japanese Island src

890

derived from seismic tomography including NIED MOWLAS Hi-net and S-net data.

891

In M. Kanao & G. Toyokuni (Eds.), Seismic Waves - Probing Earth System.

892

IntechOpen. https://doi.org/10.5772/intechopen.86936.

893

Matsumoto, K., Takanezawa, T., & Ooe, M. (2000). Ocean tide models developed by894

assimilating TOPEX/POSEIDON altimeter data into hydrodynamical model: A

895

global model and a regional model around Japan. Journal of Oceanography, 56,896

567–581. https://doi.org/10.1023/A:1011157212596

897

Miyakoshi, K., Somei, K., Yoshida, K., Kurahashi, S., Irikura, K., & Kamae, K. (2020). Scaling

898

Relationships of Source Parameters of Inland Crustal Earthquakes in Tectonically

899

Active Regions. Pure and Applied Geophysics, 177, 1917–1929.

900

https://doi.org/10.1007/s00024-019-02160-0

manuscript submitted to Journal of Geophysical Research: Solid Earth

901

Mochizuki, M., Kanazawa, T., Uehira, K., Shimbo, T., Shiomi, K., Kunugi, T., et al. (2016). S‐

902

net project: Construction of large scale seafloor observatory network for tsunamis903

and earthquakes in Japan [Conference presentation]. AGU Fall Meeting 2016,

904

American Geophysical Union, San Francisco.

905

Mulia, I. E., & Satake, K. (2021). Synthetic analysis of the efficacy of the S-net system in906

tsunami forecasting. Earth, Planets and Space, 73, 36.

907

https://doi.org/10.1186/s40623-021-01368-6

908

Murotani, S., Satake, K., & Fujii, Y. (2013). Scaling relations of seismic moment, rupture area,909

average slip, and asperity size for M~9 subduction-zone earthquakes. Geophysical910

Research Letters, 40, 5070–5074. https://doi.org/10.1002/grl.50976

911

Murotani, S. (2018) Natural disasters that reached a milestone in Fukushima –The 1888 eruption

912

of Bandai-san and the 1938 Fukushima-ken Toho-oki Earthquake– [Translated Title

913

from Japanese]. Earthquake Journal, 66, 30–39 (in Japanese). Retrieved January 16,

914

2021, from https://www.adep.or.jp/public/img/66.pdf

915

Nagai, T. (1998). Development and improvement of the Japanese coastal wave observation

916

network (NOWPHAS). Journal of Japan Society of Civil Engineers, 609, 1–14.

917

https://doi.org/10.2208/jscej.1998.609_1

918

Nakajima, J., & Hasegawa, A. (2006). Anomalous low-velocity zone and linear alignment of919

seismicity along it in the subducted Pacific slab beneath Kanto, Japan: Reactivation

920

of subducted fracture zone? Geophysical Research Letters, 33, L16309.

921

https://doi.org/10.1029/2006GL026773

922

Nakata, K., Hayashi, Y., Tsushima, H., Fujita, K., Yoshida, Y., & Katsumata, A. (2019).

923

Performance of uniform and heterogeneous slip distributions for the modeling of the924

November 2016 off Fukushima earthquake and tsunami, Japan. Earth, Planets and925

Space, 71(1), 30. https://doi.org/10.1186/s40623-019-1010-1

926
927

National Research Institute for Earth Science and Disaster Resilience [NIED] (2019). NIED S‐
net [Data set]. https://doi.org/10.17598/nied.0007

928

Nishikawa, T., Matsuzawa, T., Ohta, K., Uchida, N., Nishimura, T., & Ide, S. (2019). The slow

929

earthquake spectrum in the Japan Trench illuminated by the S-net seafloor

930

observatories. Science, 365, 808–813. https://doi.org/10.1126/science.aax5618

manuscript submitted to Journal of Geophysical Research: Solid Earth931

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

ruptures with heterogeneous stress change and implications for earthquake physics.

933

Geophysical Journal International, 193, 1691–1712.

934

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

935

Polster, A., Fabian, M., & Villinger, H. (2009). Effective resolution and drift of paroscientific936

pressure sensors derived from long-term seafloor measurements. Geochemistry,937

Geophysics, Geosystems, 10(8), Q08008. https://doi.org/10.1029/2009GC002532

938
939
940

Saito, T. (2019). Tsunami generation and propagation. Springer Japan, Tokyo.
https://doi.org/10.1007/978-4-431-56850-6
Saito, T., Satake, K., & Furumura, T. (2010). Tsunami waveform inversion including dispersive

941

waves: The 2004 earthquake off Kii Peninsula, Japan. Journal of Geophysical942

Research, 115, B06303. https://doi.org/10.1029/2009JB006884

943

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

944

Annual Review of Earth and Planetary Sciences, 48, 121–145.

945

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

946

Saito, T., Kubota, T., Chikasada, N. Y., Tanaka, Y., & Sandanbata, O. (2021). Meteorological947

tsunami generation due to sea-surface pressure change: Three-dimensional theory948

and synthetics of ocean-bottom pressure change [Preprint]. Earth and Space Science

949

Open Archive. https://doi.org/10.1002/essoar.10504961.1

950

Sasajima, R., Shibazaki, B., Iwamori, H., Nishimura, T., & Nakai, Y. (2019). Mechanism of951

subsidence of the Northeast Japan forearc during the late period of a gigantic

952

earthquake cycle. Scientific Reports, 9, 5726. https://doi.org/10.1038/s41598-019-

953

42169-y954

Sawazaki, K., & Nakamura, T. (2020). “N”-shaped Y/X coda spectral ratio observed for in-line-

955

type OBS networks; S-net and ETMC: interpretation based on natural vibration of956

pressure vessel. Earth, Planets and Space, 72, 130. https://doi.org/10.1186/s40623-

957

020-01255-6

958

Somerville, P., Irikura, K., Graves, R., Sawada, S., Wald, D., Abrahamson, N., Iwasaki, Y.,

959

Kagawa, T., Smith, N., & Kowada, A. (1999). Characterizing crustal earthquake slip

960

models for the prediction of strong ground motion. Seismological Research Letters,

961

70, 59–80. https://doi.org/10.1785/gssrl.70.1.59

manuscript submitted to Journal of Geophysical Research: Solid Earth962

Sun, T., Wang, K., Iinuma, T., Hino, R., He, J., Fujimoto, H., … Hu, Y. (2014). Prevalence of963

viscoelastic relaxation after the 2011 Tohoku-oki earthquake. Nature, 514(7520), 84–964

87. https://doi.org/10.1038/nature13778

965

Suppasri, A., Leelawat, N., Latcharote, P., Roeber, V., Yamashita, K., Hayashi, A., … Imamura,966

F. (2017). The 2016 Fukushima earthquake and tsunami: Local tsunami behavior and967

recommendations for tsunami disaster risk reduction. International Journal of

968

Disaster Risk Reduction, 21, 323–330. https://doi.org/10.1016/j.ijdrr.2016.12.016

969

Takagi, R., Uchida, N., Nakayama, T., Azuma, R., Ishigami, A., Okada, T., … Shiomi, K.

970(2019). Estimation of the orientations of the S‐net cabled ocean‐bottom sensors.

971

Seismological Research Letters, 90, 2175–2187. https://doi.org/10.1785/0220190093

972

Takagi, R., Toyokuni, G., & Chikasada, N. (2020). Ambient noise correlation analysis of S-net

973

records: extracting surface wave signals below instrument noise levels. Geophysical974

Journal International, 224, 1640–1657. https://doi.org/10.1093/gji/ggaa548

975

Tanaka, M., Asano, K., Iwata, T., & Kubo, H. (2014). Source rupture process of the 2011

976

Fukushima-ken Hamadori earthquake: How did the two subparallel faults rupture?

977

Earth, Planets and Space, 66(1), 101. https://doi.org/10.1186/1880-5981-66-101

978Tanaka, S., Matsuzawa, T., & Asano, Y. (2019). Shallow Low‐Frequency Tremor in the

979

Northern Japan Trench Subduction Zone. Geophysical Research Letters, 46, 5217–

980

5224. https://doi.org/10.1029/2019GL082817

981

Tanioka, Y. (2020). Improvement of Near-field Tsunami Forecasting Method Using Ocean

982

Bottom Pressure Sensor Network (S-net). Earth, Planets and Space, 72, 132.

983

https://doi.org/10.1186/s40623-020-01268-1

984

Terakawa, T., & Matsu’ura, M. (2010). The 3-D tectonic stress fields in and around Japan985

inverted from centroid moment tensor data of seismic events. Tectonics, 29, TC6008.

986

https://doi.org/10.1029/2009TC002626

987

Toda, S., & Goto, H. (2016)988

https://irides.tohoku.ac.jp/media/files/earthquake/eq/2016_fukushima_eq/20161122_

989

fukushima_eq_activefault_toda.pdf

990
991

Tomita, F., Kido, M., Osada, Y., Hino, R., Ohta, Y., & Iinuma, T. (2015). First measurement of
the displacement rate of the Pacific Plate near the Japan Trench after the 2011

manuscript submitted to Journal of Geophysical Research: Solid Earth992

Tohoku-Oki earthquake using GPS/acoustic technique. Geophysical Research993

Letters, 42(20), 8391–8397. https://doi.org/10.1002/2015GL065746

994

Tomita, F., Kido, M., Ohta, Y., Iinuma, T., & Hino, R. (2017). Along-trench variation in seafloor

995

displacements after the 2011 Tohoku earthquake. Science Advances, 3, e1700113.

996

https://doi.org/10.1126/sciadv.1700113

997

Tomita, F., Iinuma, T., Ohta, Y., Hino, R., Kido, M., & Uchida, N. (2020). Improvement on

998

spatial resolution of a coseismic slip distribution using postseismic geodetic data

999

through a viscoelastic inversion. Earth, Planets and Space, 72, 84.

1000

https://doi.org/10.1186/s40623-020-01207-0

1001

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

1002

inversion incorporating permanent seafloor deformation and its application to1003

tsunami forecasting. Journal of Geophysical Research, 117, B03311.

1004

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

1005

Tsushima, H., & Yamamoto, T. (2020). Operational use of tsunami source inversion in near-1006

field tsunami warning by JMA [Conference presentation]. JpGU-AGU Joint Meeting1007

2020 Virtual, Online. https://confit.atlas.jp/guide/event/jpgu2020/subject/HDS08-1008

12/advanced

1009

Uchida, N., Nakajima, J., Wang, K., Takagi, R., Yoshida, K., Nakayama, T., … Asano, Y.

1010(2020). Stagnant forearc mantle wedge inferred from mapping of shear-wave

1011

anisotropy using S-net seafloor seismometers. Nature Communications, 11.

1012

https://doi.org/10.1038/s41467-020-19541-y

1013

Uchida, N. & Bürgmann, R. (2021). A decade of lessons learned from the 2011 Tohoku-Oki

1014

earthquake. Review of Geophysics, 59, e2020RG000713.

1015

https://doi.org/10.1029/020RG000713

1016

Uehira, K., Kanazawa, T., Mochizuki, M., Fujimoto, H., Noguchi, S., Shinbo, T., et al. (2016).

1017

Outline of seafloor observation network for earthquakes and tsunamis along the1018

Japan Trench (S‐net) [Conference presentation]. EGU General Assembly 2016,

1019

European Geosciences Union, Vienna, Austria.

1020

Wallace, L. M., Araki, E., Saffer, D., Wang, X., Roesner, A., Kopf, A., … Carr, S. (2016). Near-

1021

field observations of an offshore Mw 6.0 earthquake from an integrated seafloor and1022

subseafloor monitoring network at the Nankai Trough, southwest Japan. Journal ofmanuscript submitted to Journal of Geophysical Research: Solid Earth1023

Geophysical Research: Solid Earth, 121, 8338–8351.

1024

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

1025

Wells, D. L., & Coppersmith, K. J. (1994). New empirical relationships among magnitude,1026

rupture length, rupture width, rupture area, and surface displacement. Bulletin of the1027

Seismological Society of America, 84(4), 974–1002.

1028

Wang, K., & Suyehiro, K. (1999). How does plate coupling affect crustal stresses in Northeast

1029

and Southwest Japan? Geophysical Research Letters, 26(15), 2307–2310.

1030

https://doi.org/10.1029/1999GL900528

1031

Wang, K., Sun, T., Brown, L., Hino, R., Tomita, F., Kido, M., … Fujiwara, T. (2018). Learning

1032

from crustal deformation associated with the M9 2011 Tohoku-oki earthquake.

1033

Geosphere, 14, 1–20. https://doi.org/10.1130/GES01531.1

1034

Wang, K., Brown, L., Hu, Y., Yoshida, K., He, J., & Sun, T. (2019). Stable Forearc Stressed by a1035

Weak Megathrust: Mechanical and Geodynamic Implications of Stress Changes

1036

Caused by the M = 9 Tohoku-Oki Earthquake. Journal of Geophysical Research:

1037

Solid Earth, 124, 6179–6194. https://doi.org/10.1029/2018JB017043

1038

Wang, Y., & Satake, K. (2021). Real-time tsunami data assimilation of S-net pressure gauge

1039

records during the 2016 Fukushima Earthquake. Seismological Research Letters.

1040

https://doi.org/10.1785/0220200447

1041

Watts, D. R., & Kontoyiannis, H. (1990). Deep-ocean bottom pressure measurement: drift

1042

removal and performance. Journal of Atmospheric and Oceanic Technology, 7, 296–

1043

306. https://doi.org/10.1175/1520-0426(1990)007<0296:DOBPMD>2.0.CO;2

1044

Yamamoto, N., Aoi, S., Hirata, K., Suzuki, W., Kunugi, T., & Nakamura, H. (2016a). Multi-

1045

index method using offshore ocean-bottom pressure data for real-time tsunami

1046

forecast. Earth, Planets and Space, 68, 128. https://doi.org/10.1186/s40623-016-

1047

0500-7

1048

Yamamoto, N., Hirata, K., Aoi, S., Suzuki, W., Nakamura, H., & Kunugi, T. (2016b). Rapid

1049

estimation of tsunami source centroid location using a dense offshore observation

1050

network. Geophysical Research Letters, 43, 4263–4269.

1051

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

manuscript submitted to Journal of Geophysical Research: Solid Earth1052

Yamazaki, Y., Cheung, K. F., & Lay, T. (2018). A self-consistent fault slip model for the 2011

1053

Tohoku earthquake and tsunami. Journal of Geophysical Research: Solid Earth, 123,1054

1435–1458. https://doi.org/10.1002/2017JB014749

1055

Yoshida, K., Hasegawa, A., Okada, T., Iinuma, T., Ito, Y., & Asano, Y. (2012). Stress before and1056

after the 2011 great Tohoku-oki earthquake and induced earthquakes in inland areas1057

of eastern Japan. Geophysical Research Letters, 39(3), L03302.

1058

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

1059

Yoshida, K., Hasegawa, A., & Okada, T. (2015a). Spatial variation of stress orientations in NE

1060

Japan revealed by dense seismic observations. Tectonophysics, 647, 63–72.

1061

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

1062

Yoshida, K., Hasegawa, A., & Okada, T. (2015b). Spatially heterogeneous stress field in the1063

source area of the 2011 Mw 6.6 Fukushima-Hamadori earthquake, NE Japan,

1064

probably caused by static stress change. Geophysical Journal International, 201,1065

1062–1071. https://doi.org/10.1093/gji/ggv068

1066

Yu, Z., & Zhao, D. (2020). Seismic evidence for water transportation in the forearc off northern

1067Japan. Journal of Geophysical Research: Solid Earth, 125, 2019JB018600.

1068

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

1069
1070

manuscript submitted to Journal of Geophysical Research: Solid Earth

Table 1. Fault parameters for the rectangular fault models.
Maximum vertical

Fault center location
ModelsaStrike [°] Dip [°] Rake [°]Longitude

Latitude

Depth

[°E]

[°N]

[km]a

GCMT solution

141.46

37.31

12.0

49

36

−89

Tsunami source

N/A

N/A

N/A

N/A

N/A

Grid-searchab

141.516537.3105

6.0

49

Slip distributiona

141.4908c

37.2630c

7.7c

Gusman et al. (2017)141.4532c

37.2705c

Adriano et al. (2018)141.4406c

Nakata et al. (2019)

M0 [Nm]displacement [cm]UpliftSubsidence

3.18 × 1019

N/A

N/A

N/A

N/A

16.3

238.4

36

−89

2.10 × 1019

16.0

193.1

49

36

−89

6.30 × 1019

10.5

237.4

10.1c

45

41

−95

3.70 × 1019

10.1

182.6

37.2695c

8.8c

49

35

−89

3.35 × 1019

8.5

130.6

141.4660c

37.2932c

10.5c

50

48

Variable

8.52 × 1019

29.7

222.2

JMA

141.5260c

37.2732c

10.5c

50

48

Variable

8.52× 1019

N/A

N/A

GSI

141.4971d

37.2821d

10.1d

47.6

63.2

−89.8

2.0× 1019

N/A

N/A

Fault geometry is fixed to the GCMT value.

b

Fault dimenson is L = 15 km, W = 10 km, and slip amount is D = 467.7 cm. The depths of the fault top and bottom are 3.1 km and 8.9km, respectively.
c

Slip-weighted average location is shown.

dCenter of the rectangular fault is shown. The location of the left top corner is (141.28°E, 37.17°N, 2.2 km). The fault dimension andslip amount of the rectangular fault are L = 45.1 km and W = 17.7 km, D = 0.78 m. The top and bottom depths of the rectangular fault
are 2.2 km and 18.0 km, respectively.

Journal of Geophysical Research: Solid Earth
Supporting Information for
Improving the constraint on the Mw 7.1 2016 off-Fukushima shallow normal-faulting
earthquake with the high azimuthal coverage tsunami data from the S-net wide anddense network: Implication for the stress regime in the Tohoku overriding plate
Tatsuya Kubota1, Hisahiko Kubo1, Keisuke Yoshida2, Naotaka Y. Chikasada1, Wataru Suzuki1,
Takeshi Nakamura3, and Hiroaki Tsushima41

National Research Institute for Earth Science and Disaster Resilience, Japan.
23Graduate School of Science, Tohoku University, Japan.

Central Research Institute of Electric Power Industry, Japan.

Contents of this file
Texts S1 to S3
Figures S1 to S13
Tables S1 to S3

Additional Supporting Information (Files uploaded separately)
Dataset S1 contains ocean-bottom pressure gauge data installed by Earthquake ResearchInstitute, the University of Tokyo. Dataset S2 is the tgz compressed file which contain the
results of the tsunami source modeling and the finite fault inversion.

Introduction
Texts S1, S2, and S3 explain the procedures for the tsunami source modeling, the grid-searchanalysis for the rectangular fault model, and the fault slip estimation, respectively. Figure S1 is thetrade-off curve used to determine the weight of the smoothing constraint for tsunami source
modeling. Figure S2 is the result of the tsunami source modeling using the pressure waveforms.
Figure S3 compares the tsunami source expected from the previous studies. Figure S4 is the resultof the grid-search analysis. Figure S5 evaluates the uncertainty of the fault dimension of therectangular fault model. Figure S6 is the trade-off curve used to determine the weight of thesmoothing constraint for the fault slip inversion. Figure S7 is the result of the recovery test.
Comparison with the fault models deduced from the onshore data is shown Figure S8. Figure S9shows the tsunami simulation from the previous fault models. Figures S10, S11 and S12 are theresult of the finite fault inversion only using the stations far from the focal area. Figure S13
evaluates the downdip limit of the fault depth. The time windows used for the inversion analysis issummarized in Table S1. Table S2 shows the unknown parameters searched in the grid-searchanalysis. The station locations of the OBPGs installed by Tohoku University are listed in Table S3.

1

Text S1.

This text explains the procedure for the tsunami source modeling shown in Section
4. We first explain how to simulate the tsunami Green’s function, which are the pressurechange waveforms due to the tsunami and seafloor displacement at each OBPG caused by thedisplacement of the small region of seafloor. We distribute the small elements of the seafloor
uplift (unit source elements) around the focal area (rectangular area in Figure 3a). The unit
source element of the seafloor vertical displacement is given by$$!!" (#, %) = !# (% + % cos for+%&(()(! )+$$./ 0% + % cos 1+"%&,-)-# .

+

+$
+$

#! − %" ≤ # ≤ #! + %" , %" − %$ ≤ % ≤ %" + % ,

23
(S1)

which takes the maximum value of u0 = 1 cm at (xi, yj). Here, Lx and Ly are the spatial extent of
the unit source element along the x- and y-directions, respectively. We assume that Lx = Ly = 4
km. Each of the unit sources overlaps with adjacent unit sources with a horizontal interval of
ΔLx = Lx/2 and ΔLy = Ly/2. The numbers of unit sources along the x-direction and y-directions
are Nx = 25 and Ny = 25, respectively, and the total number of unit sources is N = Nx × Ny = 625.
The size of the analytical area where the unit sources are distributed is 50 km × 50 km.
Using the seafloor vertical displacement from the unit sources, we calculate
tsunamis using the following procedure. We assume the initial sea-surface height change
assuming that the sea-surface displacement is equal to the seafloor displacement. We thensolve the linear dispersive tsunami equation (Saito et al., 2010; Saito, 2019) in Cartesian
coordinates with the staggered grid in order to simulate tsunamis:
/0

/2

$

/%

/0

/4

/2

$

/%

/0

/4

/0

/4

+ 6# ℎ /( = 3 ℎ% /(/1 - /( + /- .
/1

/4
/1

+ 6# ℎ /- = 3 ℎ% /-/1 - /( + /- .,
/2
/1

(S2)

= − /( − /-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, h is thewater depth, and g0 =9.8 m/s2 is the gravitational constant. For water depth h, we use the
JTOPO30 data with a spatial resolution of 30 arcsec, provided by the Marine Information
Research Center of the Japan Hydrographic Association (http://www.mirc.jha.jp/en/),interpolating the spatial interval of Δx = Δy = 1 km. We assume that the displacement occurs
instantaneously, at time t = 0 s. The temporal interval of the calculation is Δt = 1 s. After the
calculation, we calculate the pressure change p at each OBPG location by subtracting thepressure offset change due to the seafloor displacement from the simulated sea-surface
height change (Tsushima et al., 2012):
8 = 9# 6# (: − !5 ),(S3)where ρ0 is the density of seawater. Here we suppose seawater density ρ0 ~1.02 g/cm3 and g0=9.8 m/s2, so that a seawater column height change of 1 cm H2O can be approximated as apressure change of 1 hPa (i.e., ρ0g0 = 1 hPa/cm). We finally apply the same bandpass filter tothe simulated waveform as that applied to the observation.

2

In order to estimate the tsunami source, we use the time-derivative waveforms of
the bandpass-filtered pressure waveforms for the inversion analysis (;8/;=, Figure 3c),
because the time-derivative of the step signal becomes the impulse signal and thus does notcontain the offset change, which can reduce the artificials due to the tsunami-irrelevant steps(Kubota et al., 2018b). The data time window used for the modeling, which includes the mainpart of the tsunami (indicated by the blue traces in Figure 3c), is manually determined. We
solve the following observation equation:@>- . = - .C
AB
?

(S4)

The data vector > consists of the time-derivative waveforms of the observed pressure ;8/;=,and the matrix @ consists of the time-derivative of the tsunami Green’s functions. The vector
m consists of the amounts of the displacement of the unit sources, which are the unknown
parameters to be solved. The matrix S indicates the constraint for the spatial smoothing (e.g.,
Baba et al., 2006) and the parameter α is its weight. The goodness of the estimated source isevaluated using the variance reduction (VR):VR = F1 −∑! 78!&'( )8!)*+ 9∑! 8!&'(%%

H × 100 (%)(S5)

where M!:;< and M!=>? are the i-th data of the observed and calculated time-derivative pressure
waveforms, respectively. The smoothing weight α is determined based on the trade-off
between the weight and the VR (Figure S1) in order to avoid both the overfitting andoversmoothing of data.
Text S2.

This text explains the procedure for the grid-search analysis to estimate the planar
rectangular fault model, shown in Section 5.1. The modeling procedure is mostly based on theapproach proposed by Kubota et al. (2015; 2019). In the grid-search, we assume one planar
rectangular fault with a uniform slip. The strike angle of the fault is fixed to the GCMT value(strike = 49°), considering the consistency with the direction of the northeast-southwest
extent of the tsunami source. Since the dip and rake angles cannot be constrained only fromthe tsunami source, we assume these angles based on the GCMT solution (dip = 35° and rake =−89°), as inferred from the analysis of the teleseismic data. To find the optimum model thatbest reproduces the S-net waveforms, we vary the other fault parameters and simulate
tsunamis. The unknown parameters of the rectangular fault that we search are the fault center
location (longitude, latitude, and depth) and its dimensions (length L and width W). The slip
amount on fault D is adjusted to maximize the VR in Eq. (S5). The search range for theseparameters is summarized in Table S2, which is determined based on the tsunami source
model obtained in the previous section. Using an assumed rectangular fault with a set ofparameters (the fault model candidate), we calculated the seafloor displacement (Okada,1992). Then, using the seafloor displacemet as the initial sea surface height, we simulated
tsunamis with the same simulation scheme as the calculation of the Green's function for thetsunami source inversion (Text S1). After the calculation, we calculate the pressure changes at
the OBPG stations using Eq. (S3). Finally we evaluate the goodness of each of the fault modelcandidates is evaluated using the VR values (Eq. (S5)), using the same time window as used in
the inversion analysis (blue traces in Figure S4, Table S1).

3

Text S3.

This text describes the procedure of the finite fault slip inversion, shown in Section
5.2. The modeling procedure is almost similar to that reported by Kubota et al. (2018a). We first
assume a rectangular planar fault with dimensions of 45 km × 30 km, which is supposed topass through the optimum fault obtained by the grid search (Text S2). Then, this planar fault isdivided into small rectangular subfaults with size of 3 km × 3 km. We then simulate the Green’s
function (i.e., the pressure change waveforms excited by each subfault) using a similarcalculation procedure to that used in the grid-search analysis. The seafloor verticaldisplacements are calculated from each subfault using the equations of Okada (1992),
assuming a unit slip of 1 m. In this calculation, the strike, dip, and rake angles are fixed to the
GCMT value, as adopted in the finite fault inversion. Then tsunami is calculated using thevertical displacement distribution as the initial sea surface height change (Eq. (S2)). After the
calculation, the pressure change is calculated by subtracting the component of the seafloor
vertical deformation from the simulated sea surface height change (Eq. (S3)), to obtain theGreen's function for the finite fault slip inversion.
Using the Green's function for the finite fault slip inversion, simulated by theprocedure shown above, we solve the following observation equation, which is similar to thetsunami source inversion:N>- . = - . C.
AB
?

(S6)

The data vector > consists of the time-derivative waveforms of the observed pressure ;8/;=,and the matrix N consists of the time-derivative of the Green's functions for the finite fault slip
inversion. The vector m consists of the fault slip amount for each subfault (unit: [m]), which are
the unknown parameters to be solved. The matrix S indicates the constraint for the spatialsmoothing (e.g., Baba et al., 2006) and the parameter α is its weight.
When solving this observation equation, we imposed a nonnegativity constraint(Lawson & Hanson, 1974) because the negative slip (i.e., reverse-faulting slip component) is
quite unlikely to occur. The weighting of the smoothing constraint α is determined based on
the trade-off curve between its weight and VR value (Figure S6).

4

Figure S1. Trade-off curve between the smoothing weight α and VR. Red and blue solid lines
are the trade-off curves for the inversions using the time-derivative waveform of the pressure(Figure 4) and the pressure waveform (Figure S2), respectively. Dashed lines denote theweight values used for the inversion analyses. Note that the VR values are calculated by using
the time-derivative pressure waveforms for the inversion using the time derivative-waveforms(red), and by using the pressure waveforms for the inversion using the pressure waveforms(blue), respectively.

5

(b) Pressure(c) Time-derivative of pressure

S3N20
S3N21
S3N22
S3N24
S3N25
S3N26

(a)
GCMT

S2N01
S2N02
S2N03
S2N04
S2N05
S2N06
S2N07
S2N08
S2N09
S2N10

37.5°

S2N11
S2N12
S2N13

50 hPa1 hPa/s

S2N14
S2N15

37°
0

10 20
[cm]

Contour:20cm
141°−200−100 0
141.5°

100 200
142°

S2N16
S2N17
S2N18
S2N19
S2N20
S2N21
S2N22
S2N23
S2N24
S2N25
S2N26
S1N01
S1N02
S1N03
S1N04
S1N05

−10 0 10 20 30 40 50 60 −10 0 10 20 30 40 50 60

Time [min]Time [min]Figure S2. Results of the tsunami source inversion based on the conventional method. (a) Spatial distribution of the tsunami source. Contour
lines denote the subsided region with intervals of 20 cm. The white star is the JMA epicenter. The yellow and gray circles show the S-net
OBPGs used or not used, respectively, for inversion analysis. Comparisons of (b) the pressure waveforms and (c) the time-derivative
waveforms. The gray and red traces denote the observed waveforms and simulated waveforms from the tsunami source model. Traces
marked by blue lines denote the time window used for the inversion analysis.

1(a) Gusman et al. (2017)(b) Adriano et al. (2018)(c) Nakata et al. (2019)

37.5°

37°
0

10 20

0

10 20

[cm]
−200−100 0
141°141.5°

100 200
142°

0

10 20

[cm]
−200−100 0
141°141.5°[cm]

100 200
142°−200−100 0
141°141.5°

100 200
142°

Figure S3. Comparison of the tsunami source calculated from the finite fault models of the previous studies (black contours) and the tsunami
source model (red). Models of (a) Gusman et al. (2017), (b) Adriano et al. (2018), and (c) Nakata et al. (2019) are shown. The contour intervals
are 20 cm. The configuration of the fault is also shown by gray lines.

2(b) Pressure(c) Time-derivative of pressure

S3N20
S3N21
S3N22
S3N23
S3N24
S3N25
S3N26

(a)
GCMT

S2N01
S2N02
S2N03
S2N04
S2N05
S2N06
S2N07
S2N08
S2N09
S2N10

37.5°

S2N11
S2N12
S2N13

50 hPa1 hPa/s

S2N14
S2N15

37°
0

10 20
[cm]

Contour:20cm
141°−200−100 0
141.5°

100 200
142°

S2N16
S2N17
S2N18
S2N19
S2N20
S2N21
S2N22
S2N23
S2N24
S2N25
S2N26
S1N01
S1N02
S1N03
S1N04
S1N05

−10 0 10 20 30 40 50 60 −10 0 10 20 30 40 50 60

Time [min]Time [min]Figure S4. Results of the grid-search analysis. (a) Spatial distribution of the tsunami source. The green rectangle shows the location of therectangular fault model. Black contours are the tsunami source calculated from the rectangular fault model (Table 1). The distribution of thetsunami source model obtained by the inversion is also shown by gray contours. Comparisons of (b) the pressure waveforms and (c) the timederivative waveforms. See Figure 3 and S2 for a detailed explanation of the figure.

3

141°

141.5°

142°

(a)
GCMT

37.5°

37°
0

10 20
[cm]

−200−100 0

Contour:20cm(b) S2N01

100 200(c) S2N02(d) S2N09(e) S2N10

W=10km

L=40kmL=35km
L=30km
L=25km
L=20km
L=15km
L=10km
L=5km
30 hPaL=15km20 hPa10 hPa10 hPaW=30km
W=25km
W=20km
W=15km
W=10km
W=5km(f) S2N11(g) S2N12(h) S2N14(i) S2N15

W=10km

L=40kmL=35km
L=30km
L=25km
L=20km
L=15km
L=10km
L=5km
20 hPaL=15km30 hPa40 hPa40 hPaW=30km
W=25km
W=20km
W=15km
W=10km
W=5km−10 0 10 20 30 40 50 60 −10 0 10 20 30 40 50 60 −10 0 10 20 30 40 50 60 −10 0 10 20 30 40 50 60

Time [min]Time [min]Time [min]Time [min]Figure S5. (a) Horizontal location of the optimum rectangular fault model. The green
rectangle shows the location of the rectangular fault model and countour lines are the
distribution of the tsunami source calculated from the rectangular fault model. (b–i)
Evaluation of the fault dimension. Comparisons of the stations near the epicenter between theobserved (gray) waveform and the simulated waveforms from the varied fault dimensions areshown. The simulated waveforms with thick red and blue traces denote the optimum
rectangular fault.

1Figure S6. Trade-off curve between the smoothing weight α and VR for the finite fault slip
inversion. Red solid line is the trade-off curves for the inversions using the time-derivative
waveform of the pressure (Figure 4). Dashed lines denote the weight values adopted for the
inversion analyses.

2

Max Slip = 18.2cm(a) 3×3 km

Max Slip = 46.9cm(b) 6×6 km

Max Slip = 101.3cm(c) 9×9 km

37.5°

37°
0
0

10 20

0
[cm]
100

50

Max Slip = 112.4cm(d) 12×12 km

0

10 20

0
[cm]
100

50

Max Slip = 120.3cm(e) 15×15 km

0

10 20
[cm]
100

50

Max Slip = 117.0cm(f) 15×12 km

37.5°

37°
0
0
141°

141.5°

10 20
50

0
[cm]
100

142°

0
141°

141.5°

10 20
50

0
[cm]
100

142°

0
141°

141.5°

10 20
50

[cm]
100
142°

Figure S7. Results of the recovery test of the inversion analysis. The assumed faults sizes are (a) 3 × 3 km, (b) 6 × 6 km, (c) 9 × 9 km, (d) 12 × 12km, (d) 15 × 15 km, and (e) 15 × 12 km, respectively. The recovered slip distribution is shown by colors, and black rectangles denote thelocation of the assume fault.

1(a) JMA teleseismic model (Maximum slip = 4.02 m)(b) GSI onshore geodetic model (Slip=0.78m)

37.5°

37°
0

10 20

0

10 20

Contour:1m
141°

141.5°

142°

141°

141.5°

142°

Figure S8. Comparison with the fault model estimated from the onshore data. (a) Comparison with the fault model obtained by using theteleseismic data provided by JMA (https://www.data.jma.go.jp/svd/eqev/data/sourceprocess/event/2016112205594689far.pdf). The gray andpink contour lines denote the slip distribution of the 2016 off-Fukushima earthquake obtained by this study and JMA, respectively (contour
interval: 1 m). Large dark gray rectangles show the configuration of the fault plane of the JMA analysis.
(b) Comparison with the fault model obtained by using the onshore geodetic data provided by Geospatial Information Authority of Japan(GSI, https://cais.gsi.go.jp/YOCHIREN/activity/214/214.e.html, https://cais.gsi.go.jp/YOCHIREN/activity/214/image214/008.pdf). The blue
rectangle denotes the location of the rectangular planar fault model estimate by GSI. The parameters are follows: fault left top corner location= (37.17°N, 141.28°E, 2.2km), length = 45.1 km, width = 17.7 km, strike = 47.6°, dip = 63.2°, rake = –89.8°, and slip amount = 0.78 m.

2

(a) North stations (~38°N)
S3N20

S3N24

S3N25

S2N01

Inversion
Gusman et al. (2017)
Adriano et al. (2018)
Nakata et al. (2019)
5 hPa

20

30

10 hPa40

20Time [min]30

15 hPa40

20Time [min]3030 hPa40

20Time [min]30

40Time [min](b) Middle stations (~36.5-37°N)
S2N10

S2N15

S2N16

S2N17

Inversion
Gusman et al. (2017)
Adriano et al. (2018)
Nakata et al. (2019)15 hPa2030

45 hPa40

20Time [min]3040 hPa40

20Time [min]3030 hPa40

20Time [min]30

40Time [min](c) Southern stations (~36°N)
S1N01

S2N26

S1N02

S1N03

Inversion
Gusman et al. (2017)
Adriano et al. (2018)
Nakata et al. (2019)20 hPa20

30Time [min]25 hPa40

20

30Time [min]15 hPa40

20

30Time [min]10 hPa

40

20

30

40

Time [min]Figure S9. Comparison of the observed (gray) and synthesized S-net OBPG waveforms, fromthe fault models estimated by this study (red) and the previous studies (blue, Gusman et al.,2017; Adriano et al., 2018; Nakata et al., 2019). (a) Comparison for the stations at northern partof the off northeastern Japan (~38 °N). (b) Comparison for the stations at middle part of the off
northeastern Japan (~36.5–37 °N). (c) Comparison for the stations at middle part of the off
northeastern Japan (~36 °N).

1

(a) Slip distribution
GCMT(b) Shear stress change
GCMT

37.5°

37°
0
0
141°et al. (2017)(c) Gusman
GCMT

10 20
350

141.5°

0
[cm]
700

142°

−10 −5141°et al. (2018)(d) Adriano
GCMT

10 20
0

141.5°

[MPa]
10

5

142°

(e) Nakata et al. (2019)
GCMT

37.5°

37°
0
0
141°

141.5°

10 20
350

0
[cm]
700

142°

0
141°

141.5°

10 20
350

0
[cm]
700

142°

0
141°

141.5°

10 20
350

[cm]
700
142°

Figure S10. Result of the slip inversion using only the stations far from the focal area. (a) Slip distribution (colored tiles). The pink and whitestars indicate the slip-weighted averaged centroid and the JMA epicenter, respectively. (b) Shear stress change along the fault. Negative(blue) and positive (red) denote the shear stress decrease and increase, respectively. The dark red rectangle denotes the optimum rectangularfault obtained by the grid-search analysis. (c–e) Slip distribution by the previous studies.

1(b) S-net OBPG(a)(c) Time-derivative

S3N20
S3N21
S3N22
S3N23
S3N24
S3N25
S3N26
S2N01
S2N02
S2N03
S2N04
S2N05
S2N06
S2N07
S2N08
S2N09
S2N10
S2N11
S2N12

GCMT

37.5°

S2N13

50 hPa1 hPa/s

S2N14

37°
0

10 20
[cm]

Contour:20cm
141°−200−100 0
141.5°

100 200
142°

S2N15
S2N16
S2N17
S2N18
S2N19
S2N20
S2N21
S2N22
S2N23
S2N24
S2N25
S2N26
S1N01
S1N02
S1N03
S1N04
S1N05

−10 0 10 20 30 40 50 60 −10 0 10 20 30 40 50 60

Time [min]Time [min]Figure S11. Result of the slip inversion using only the stations far from the focal area. (a) Slip distribution (colored tiles) and tsunami source
distribution (black contours, 20 cm interval). (b–c) Comparisons of the observed (gray) and simulated (red) waveforms.

2

Max Slip = 7.3cm(a) 3×3 km

Max Slip = 23.0cm(b) 6×6 km

Max Slip = 51.2cm(c) 9×9 km

37.5°

37°
0
0

10 20

0
[cm]
100

50

Max Slip = 70.6cm(d) 12×12 km

0

10 20

0
[cm]
100

50

Max Slip = 84.0cm(e) 15×15 km

0

10 20
[cm]
100

50

Max Slip = 76.3cm(f) 15×12 km

37.5°

37°
0
0
141°

141.5°

10 20
50

0
[cm]
100

142°

0
141°

141.5°

10 20
50

0
[cm]
100

142°

0
141°

141.5°

10 20
50

[cm]
100
142°

Figure S12. Result of the recovery test, which used only the stations far from the focal area (Figure S10). The other caption is identical to
Figure S7.

3

141°

141.5°

142°

(a)
37.5°
km
.2 m
11 4.6k km
1 8.0
1

37°
0

10 20

(b) North stations (~38°N)
Inversion
Zbot=9.4km
Zbot=11.2km
Zbot=12.9km
Zbot=14.6km
Zbot=16.3km
Zbot=18.0km

S3N25

S3N26

S2N01

S2N02

15 hPa20 hPa30 hPa20 hPa−10 0 10 20 30 40 50 60 −10 0 10 20 30 40 50 60 −10 0 10 20 30 40 50 60 −10 0 10 20 30 40 50 60

Time [min]Time [min]Time [min]Time [min](c) Middle stations (~36.5°N)
Inversion
Zbot=9.4km
Zbot=11.2km
Zbot=12.9km
Zbot=14.6km
Zbot=16.3km
Zbot=18.0km

S2N16

S2N17

S2N18

S2N19

30 hPa15 hPa10 hPa10 hPa−10 0 10 20 30 40 50 60 −10 0 10 20 30 40 50 60 −10 0 10 20 30 40 50 60 −10 0 10 20 30 40 50 60

Time [min]Time [min]Time [min]Time [min](d) Southern stations (~36°N)
Inversion
Zbot=9.4km
Zbot=11.2km
Zbot=12.9km
Zbot=14.6km
Zbot=16.3km
Zbot=18.0km

S2N26

S1N01

S1N02

S1N03

20 hPa25 hPa15 hPa10 hPa−10 0 10 20 30 40 50 60 −10 0 10 20 30 40 50 60 −10 0 10 20 30 40 50 60 −10 0 10 20 30 40 50 60

Time [min]Time [min]Time [min]Time [min]Figure S13. The simulation of the simple rectangular fault with different fault widths. (a)Configurations of the assumed faults. Depths of the bottom end of the rectangular faults (Zbot)
are 9.4, 11.2, 12.9, 14.6, 16.3, and 18.0 km. (b–d) Comparisons of the observed (gray)waveforms and the waveforms simulated from the finite fault model (red) and simulated fromthe assumed rectangular fault models (blue).

1

Table S1. Time window range used for the inversion analysis.
Begin time of time
End time of time
Stationdt [s]
awindow [s]window [s]a
S3N20
1200
2700
10
S3N21
1200
2700
10
S3N22
900
2400
10
S3N23
900
2400
10
S3N24
900
2400
10
S3N25
900
2400
10
S3N26
1500
3000
10
S2N01
900
2400
10
S2N02
300
1800
10
S2N03
600
2100
10
S2N04
600
2100
10
S2N05
900
2400
10
S2N06
900
2400
10
S2N07
900
2400
10
S2N08
600
2100
10
S2N09
600
2100
10
S2N10
600
2100
10
S2N11
300
1800
10
S2N12
150
1650
10
S2N13
N/A
N/A
N/A
S2N14
150
1650
10
S2N15
600
2100
10
S2N16
600
2100
10
S2N17
600
2100
10
S2N18
600
2100
10
S2N19
600
2100
10
S2N20
600
2100
10
S2N21
900
2400
10
S2N22
900
2400
10
S2N23
900
2400
10
S2N24
900
2400
10
S2N25
900
2400
10
S2N26
900
2400
10
S1N01
1500
3000
10
S1N02
1200
2700
10
S1N03
1200
2700
10
S1N04
1200
2700
10
S1N05
1200
2700
10
a
Time window is measure from the origin time.
b
In the inversion, the bandpass-filtered (100–3600 s) data are resampled every dt s whensolving the observation equation.

2

Table S2. Search range for the grid search analysis.
Parameters
Range
IncrementabLongitude141.46°E ± 20 km
5 km
Latitude ab
37.31°N ± 20 km
5 km
Depth ab12.0 km ± 10 kma
2 kma
Strike
49°
Fixed
Dip a
35°
Fixed
Rake a−89°
Fixedc
Length
5 km – 60 km
5 km
Widthc5 km – 60 km
5 km
Slip amount
Adjusted so that the VR value takes the maximumaReference values are taken from the GCMT solution.
b
Fault center location is shown.
c
When the depth of the updip end of the fault is shallower than a depth of 0.1 km, thecalculation is skipped.
Table S3. Station list of the OBPGs installed by Tohoku University
Observation duration
LoggerStationLongitude (°E)
Latitude (°N)
Depth (m)(yyyy/mm/dd)
typea
G01
144.920438.70305456
2016/05/22 – 2017/04/11
UME
G12
143.5317
38.02134366
2016/05/24 – 2017/04/10
UME
G16
143.0470
37.33244414
2016/05/27 – 2017/04/15
HAK
G17b
142.7123
36.89794232
2016/05/28 – 2017/04/09
HAK
G19
142.6735
36.49315691
2016/05/28 – 2017/04/09
HAK
AoA50
142.3176
36.87252853
2016/09/22 – 2017/11/09
UME
AoA60b142.7140
36.89934225
2016/09/22 – 2017/10/15
UME
AoA70
142.2868
36.69372544
2016/09/22 – 2017/10/15
HAK
FKOB1
142.5800
36.80554550
2016/09/28 – 2017/10/15
UME
FKOB2
142.8553
36.72255506
2016/09/28 – 2017/10/14
HAK
G15
143.5215
37.67735239
2016/10/02 – 2017/10/19
UME
a
UME: Paroscientific Series 8CB intelligent type pressure sensor + Umezawa-Musen Co. datalogger, HAK: Paroscientific Series 8B pressure sensor + Hakusan Co. LS9150 data loggerbStation G17 and AoA60 are installed at almost identical location.

3