2.1 Modeling Strategy
To understand the mechanics of the shallow trench-reaching rupture, it is essential to reveal the fault rupture mechanics (more specifically, the coseismic stress drop at the shallow portion) and to precisely estimate the shallow slip with ultra-high resolution. Seafloor geodetic data (Fujiwara et al., 2011; Kido et al., 2011; Sato et al., 2011) have made large contributions towards elucidating this issue (Iinuma et al., 2012; Sun et al., 2017); however, the spatial resolution of such data is limited, which contains information that is only relevant at the location of the observation point. In addition to these offshore, and onshore, geodetic datasets (Fig. 1c), the present study utilized tsunami waveform data (Figs. 1a and 1b) which contain unique and robust spatial information about the distribution of the fault slip, in particular at the shallow portion near the trench, as the spatial extent of the tsunami source. This is because shallow near-trench slips excite tsunamis more efficiently than the deeper slip. Furthermore, we used data from ocean-bottom pressure gauges installed just above the large slip area (Kubota et al. 2021), which have not been used before in any modeling analyses (Fig. 1b). The advantages of combining the far-field tsunami waveforms and the in-situ ocean-bottom pressure gauge waveforms contribute to a significant improvement of the fault model construction, particularly the estimation with the high-resolution on the fault slip at the shallow portion. In contrast, the onshore datasets, which have been widely used in past studies but obtained far from the focal area had a worse constraint on the shallow slip. The onshore geodetic data, recorded more than 200 km from the trench axis, have difficulty in resolving the slip near the trench (Loveless and Meade, 2015). The analyses using the onshore seismograms require the assumption of the rupture propagation velocity across the fault, but this assumption may cause a large uncertainty in estimating the spatial extent of the entire slip region because of a substantial trade-off between the rupture velocity and the spatial extent (Kubota et al., 2018).
Some conventional fault modeling methods using tsunami data (e.g., Satake et al., 2013; Yamazaki et al., 2018) have often simplified the fault geometry by assuming planar subfaults, which do not reflect the real three-dimensional fault geometry. This simplification may cause serious problems in inferring the correct fault slip and stress changes on the fault (Wang et al., 2018). Thus, we incorporated the non-planar plate geometry model (Koketsu et al., 2012) using small triangular subfault elements. Note that the shallowest part of the fault geometry was slightly modified from the plate geometry model, to appropriately simulate trench-breaching rupture (Wang et al., 2018, see Section 2.3.1 for details).
2.2 Data
We used the tsunami data recorded by the ocean-bottom pressure gauges installed around the focal area, obtained by Tohoku University (Kubota et al., 2021) (green inverted triangles in Figs. 1a and 1b) and the University of Tokyo (Maeda et al., 2011) (dark red inverted triangles in Fig. 1b). The station information is summarized in Table 1. After removing the signal fluctuation due to the ocean tide, we applied a low-pass filter with a cut-off period of 100 s to these records.
Table 1
List of tsunami stations used in this study
Station | Latitude [°N] | Longitude [°E] | Depth [m] | Inversion time window [s] | Agency | Sampling rate of original data [s]a |
TM2 | 39.2459 | 142.4526 | 997 | 0–1800 | ERI | 0.1 |
TM1 | 39.2283 | 142.7720 | 1618 | 0–1800 | ERI | 0.1 |
P06 | 38.6340 | 142.5838 | 1254 | 0–3600 | Tohoku University | 1 |
P02 | 38.5002 | 142.5016 | 1104 | 0–3600 | Tohoku University | 1 |
P03 | 38.1834 | 142.3998 | 1052 | 0–3600 | Tohoku University | 1 |
P07 | 38.0016 | 142.4495 | 1059 | 0–3600 | Tohoku University | 1 |
P08 | 38.2829 | 142.8320 | 1418 | 0–3600 | Tohoku University | 1 |
P09 | 38.2650 | 143.0002 | 1556 | 0–3600 | Tohoku University | 1 |
GJT3 | 38.2945 | 143.4814 | 3293 | 0–3600 | Tohoku University | 1 |
21418 | 38.7180 | 148.6980 | 5500 | 1200–4200 | DART | 15 |
KPG2 | 42.2365 | 144.8454 | 2210 | 600–3600 | JAMSTEC | 0.1 |
KPG1 | 41.7040 | 144.4375 | 2218 | 600–3600 | JAMSTEC | 0.1 |
KCTD | 41.6675 | 144.3409 | 2540 | 600–3600 | JAMSTEC | 10 |
NMS09 | 42.3692 | 145.9167 | 3316 | 600–3600 | Tohoku University | 1 |
NMS05 | 42.1667 | 145.8235 | 4548 | 600–3600 | Tohoku University | 1 |
BOSO2 | 34.7550 | 140.7517 | 2098 | 600–4200 | JMA | 1 |
BOSO3 | 34.8050 | 140.5067 | 1912 | 600–4200 | JMA | 1 |
HPG1 | 35.0031 | 139.2247 | 1176 | 1800–5400 | JAMSTEC | 1 |
VCM3 | 35.0712 | 139.3906 | 1225 | 1800–5400 | NIED | 0.1 |
VCM1 | 34.5954 | 139.9198 | 2125 | 1800–5400 | NIED | 0.1 |
807 | 40.1167 | 142.0667 | 125 | 0–3600 | NOWPHAS | 5 |
804 | 39.6272 | 142.1867 | 200 | 0–3600 | NOWPHAS | 5 |
802 | 39.2586 | 142.0969 | 204 | 0–3600 | NOWPHAS | 5 |
803 | 38.8578 | 141.8944 | 160 | 0–3600 | NOWPHAS | 5 |
801 | 38.2325 | 141.6836 | 144 | 0–3600 | NOWPHAS | 5 |
806 | 36.9714 | 141.1856 | 137 | 0–3600 | NOWPHAS | 5 |
613 | 42.9106 | 144.3972 | 50.0 | 1800–5400 | NOWPHAS | 5 |
602 | 42.5439 | 141.4458 | 50.7 | 1800–5400 | NOWPHAS | 5 |
202 | 40.9250 | 141.4242 | 43.8 | 1800–4800 | NOWPHAS | 5 |
203 | 40.5608 | 141.5683 | 27.7 | Not Used | NOWPHAS | 5 |
219 | 40.2178 | 141.8600 | 49.5 | Not Used | NOWPHAS | 5 |
205 | 38.2500 | 141.0661 | 21.3 | Not Used | NOWPHAS | 5 |
aAll observed records were resampled to 1 s in the inversion analyses. |
We also used the offshore tsunami data obtained far from the source area (far-field data) from the ocean-bottom pressure gauges of the Japan Agency for Marine-Earth Science and Technology (JAMSTEC, light blue inverted triangles in Fig. 1a), Japan Meteorological Agency (JMA, red inverted triangles), National Research Institute for Earth Science and Disaster Resilience (NIED, black inverted triangles), and National Oceanic and Atmospheric Administration (NOAA)’s DART (Deep-ocean Assessment and Reporting of Tsunamis) system (a blue inverted triangle). We also used the tsunami waveforms from the GPS buoys (yellow squares in Figs. 1a and 1b) and the wave gauges (orange triangles) of NOWPHAS (Nationwide Ocean Wave information network for Ports and HArbourS) of the Port and Airport Research Institute (PARI). The data processing was similar to the pressure data around the focal area, but a bandpass filter with a passband of 100–3600 s was used. We also used onshore GPS data obtained by the Geospatial Information Authority of Japan (GSI), as well as offshore geodetic observation data (Fig. 1c) (Kido et al., 2011; Sato et al., 2011).
2.3 Methods
2.3.1 Fault Geometry and Crustal Deformation Calculation
As mentioned above, conventional fault modeling using tsunami data has often simplified the fault geometry by using methods such as implementing large-sized planar faults and/or a buried fault whose top does not reach the free surface. However, these simplifications, which do not reflect the real three-dimensional fault geometry, may cause serious problems when inferring fault slip (Wang et al., 2018). Serious problems also occur in the calculation of the stress drop distribution. If large-sized subfaults with uniform slip are used, the stress concentration and/or the artificial discontinuity of the stress drop occur at the boundaries of the subfaults. If the fault top does not reach the free surface, the spurious stress increase occurs in the area between the free surface and the buried fault top. Hence, it is essential to appropriately consider the three-dimensional fault geometry to obtain the correct stress drop distribution.
In contrast to most of the past tsunami modeling studies using the simple planar fault configuration, we incorporated the nonplanar fault geometry of the Japan Integrated Velocity Structure Model (JIVSM) (Koketsu et al., 2012) for the analysis. To derive the fault slip distribution via Green's functions, the plate geometry was divided into small triangular subfault elements, in which the length of one side of the triangle was approximately 10 km, referring to the subfault configuration of a previous study (Iinuma et al., 2016) (triangles in Fig. 1). Using these triangular subfault elements, we calculated the seafloor displacement using a uniform half-space structure model (Meade, 2007). To correctly simulate the trench-extending rupture (Wang et al., 2018), we slightly modified the fault configuration as follows: First, the depths of all the triangular vertices were systematically moved vertically towards the surface by 8 km, considering the seawater depth around the trench axis. Then, the depths of the uppermost row of triangular subfaults (i.e. the computational free surface) were set to 0 km so that the shallowest fault surface accorded with the trench axis.
2.3.2 Estimation of Fault Slip Distribution
The procedure for the analyses to estimate the slip distributions (Fig. 2a) was similar to that used in our previous studies (Kubota et al., 2021). Considering that the observed data is expressed by the linear superposition of Green's function, we estimated the amount of slip at each of the triangular subfaults by solving the following equation:
$$\left(\begin{array}{c}{w}_{\text{t}\text{s}\text{u}\text{n}}{\mathbf{d}}_{\text{t}\text{s}\text{u}\text{n}}\\ {w}_{\text{g}\text{e}\text{o}\text{d}}{\mathbf{d}}_{\text{g}\text{e}\text{o}\text{d}}\\ 0\\ 0\end{array}\right)=\left(\begin{array}{c}{w}_{\text{t}\text{s}\text{u}\text{n}}{\mathbf{G}}_{\text{t}\text{s}\text{u}\text{n}}\\ {w}_{\text{g}\text{e}\text{o}\text{d}}{\mathbf{G}}_{\text{g}\text{e}\text{o}\text{d}}\\ \alpha S\\ \beta E\end{array}\right)\mathbf{m}$$
1
.
The vector d is the data vector, consisting of the observed data, and G is the matrix consisting of Green’s functions (the synthetics by the unit slip of the subfault). The subscripts denote the types of datasets. The scalar value w denotes the weight of each dataset. The weight of the tsunami data is ten times larger than the geodetic data. Vector m consists of the slip amount of the triangular subfaults, which is to be estimated. To stabilize the solution, we imposed the smoothing constraint (Maerten et al., 2005) expressed by matrix S and the damping constraint using the identity matrix E. Parameters α and β are the weights of each constraint.
The procedure to calculate Green’s functions for tsunamis (matrix Gtsun in Eq. (1)) is as follows: The initial seafloor vertical displacement distribution from each of the triangular subfaults (Meade, 2007) was simulated. The total number of the triangular subfaults used in this study was Nsub = 434. The effect of the apparent seafloor vertical movement due to the horizontally moving seafloor was also incorporated (Tanioka and Satake, 1996). We then simulated the initial sea-surface height change distribution from the seafloor deformation using a spatial smoothing filter related to the seawater depth (Saito, 2019). Finally, we simulated tsunamis using the linear dispersive equation (Saito, 2019). The cosine-shaped source time function with the duration of Tr = 40 s is assumed for the rupture time history of each subfault, as used in Kubota et al. (2021). In this simulation, the bathymetry data of GEBCO 2020 were interpolated to 2 km spatial intervals. The time interval for this simulation was 1 s. Green’s functions for the geodetic data (Ggeod) were also calculated similarly from each of the triangular subfaults (Meade, 2007).
We consider the temporal evolution of the rupture by distributing the Green’s function in the time domain in addition to the space domain. We distributed the Green’s functions in the time domain for each subfault with the temporal interval of Δt = 20 s (i.e., the slip of the k-th temporal element begins at t = tkbeg = (k − 1) × Δt and ends at t = tkend = tkbeg+ Tr). We assumed Nt = 9 Green’s functions in the time domain for each subfault. The total number of unknown parameters was N = Nsub × Nt = 3906. Considering the rupture front propagation, the slips of the k-th temporal element at the i-th subfault were constrained to be zero when the rupture front did not arrive there (i.e., we allowed the slip only when satisfying the condition Vr × tkend ≥ ai, ai is the distance between the hypocenter and the center of the i-th subfault and Vr = 4 km/s is the rupture velocity).
Using Green’s function, Eq. (1) was solved. The time windows used for this inversion analysis are shown in the blue traces in Fig. 3a–d. The weights of the spatial smoothing and damping constraints (parameters α and β) were set based on the previous modeling results (Kubota et al., 2021).
2.3.3 Evaluation of Fault Slip
We evaluated our fault slip distribution model (Fig. 2a) based on several approaches. First, the modeling uncertainty of the shallow slip estimation was evaluated by an inversion analysis based on the jack-knife or leave-one-out approach (red shaded area in Fig. 4a, in the Off-Miyagi region). In this test, we excluded one of the nine near-field pressure gauge stations and used the remaining eight stations to estimate the slip distribution. The inversion setting was identical to the original setting, which used all the data. The possible range of the slip amount was defined by the maximum and minimum slip amounts from all tested models.
We also examine the improvement in the slip distribution of our model by conducting tsunami simulations using the previous slip distributions of Satake et al. (2013), Iinuma et al. (2012), and Yamazaki et al. (2018) (Fig. 5). For the model of Satake et al. (2013) and Yamazaki et al. (2018), the propagation of rupture front was considered as done in these studies. The multiple source functions in the time domain were taken into account in the model of Satake et al. (2013) and the single source time function with duration of Tr = 32 s in the time domain was assumed in the model of Yamazaki et al. (2018). The temporal evolution was not considered in the static fault model of Iinuma et al. (2012), but we assumed that the slips of all subfaults begin simultaneously at t = 0 with the rise time of Tr = 60 s. The tsunami simulations were conducted using the linear long equation (Saito, 2019).
We further examined whether the observed near-field pressure waveforms could be explained by the slip profile estimated from the change in the bathymetry near the trench (Sun et al., 2017) (Fig. 6). We assumed the slip amounts of the i-th triangular subfaults (\({m}_{i}^{\text{'}}\)) to emulate the along-dip slip profile in Sun et al. (2017), based on the following formula:
$${m}_{i}^{\text{'}}={m}_{i}+{\Delta }{m}_{i}$$
3
.
Here, \({m}_{i}\) is the slip amount of the present model (Fig. 1b) and \({\Delta }{m}_{i}\) is the amount of modification. The modification value is defined as:
$${\Delta }{m}_{i}=\left\{\begin{array}{c}\begin{array}{cc}{d}_{0}& \left(0\le {r}_{i}<{r}_{1}\right)\end{array}\\ \begin{array}{cc}2{d}_{0}-{d}_{0}\frac{{r}_{i}}{{r}_{2}-{r}_{1}}& \left({r}_{1}\le {r}_{i}<{r}_{2}\right)\end{array}\\ \begin{array}{cc}0& \left({r}_{i}\le {r}_{2}\right)\end{array}\end{array}\right.$$
4
,
where ri is the horizontal distance between the center location of the i-th triangular subfault, xi = (xi, yi), and the reference point, x0 = (x0, y0) = (144.0°E, 38.08°N), defined as:
$${r}_{i}=\sqrt{{\left({x}_{i}-{x}_{0}\right)}^{2}+{\left({y}_{i}-{y}_{0}\right)}^{2}}$$
5
We used the values of r1 = 50 km, r2 = 100 km, and d0 = 12.5 m. The modification value with a function of the horizontal distance from point x0 is shown by a grey line in Fig. 6a.
2.3.4 Calculation of Stress Drop Distribution
After estimating the slip distribution of the triangular subfaults, we calculated the distribution of the shear stress change along the fault (i.e. stress drop, Fig. 2b) by computing the shear stress change along the slip direction at the center of each subfault:
$${\Delta }{\sigma }_{i}={\sum }_{i}{\Delta }{\sigma }_{ij}^{0}{m}_{j}$$
2
,
where \({\Delta }{\sigma }_{i}\) is the stress drop at the i-th fault, \({\Delta }{\sigma }_{ij}^{0}\) is the stress drop at the center of the i-th triangular subfault by the unit slip at the j-th subfault, and mj is the slip amount at the j-th subfault.
2.3.5 Evaluation of stress drop area
To validate the location of the large stress drop area and to examine the stress drop amount at the shallow portion, we conducted additional tsunami simulations (Figs. 7 and 8). First, we constructed the fault slip distribution models of the triangular subfaults. In this procedure, we assigned a stress drop of 5 MPa in the shallow, deep, and deeper portions, respectively (Figs. 7a–7c, shown by thick black lines). Then, considering the given stress drop amount as the data (left-hand side of Eq. (1)), the slip amount of the j-th subfault (\({m}_{j}\)) was estimated by solving the linear inversion problem, and then tsunamis were calculated (Fig. 7d).