We conducted numerical simulations of earthquake generation cycles using the realistic three-dimensional geometry of the subducting Pacific Plate along the Japan Trench30. We used the same equations, initial conditions, seismic radiation damping term, plate geometry, and plate convergence rate as in our previous study31, which roughly approximated several characteristics of M7 earthquakes (such as the 1978 Miyagi–Oki earthquake) and M9 earthquakes (such as the 2011 Tohoku-Oki earthquake) in the middle segment of the Japan Trench. Seismic and aseismic events were modeled to represent the release of the slip deficit or backslip that accumulates during interseismic periods32. Such spatiotemporal variations in the slip velocity were assumed to indicate an unstable slip with a frictional interface. We used a rate- and state-dependent frictional law16 as an approximate mathematical model for large-scale frictional behavior on the plate interface and a fault constitutive law33 to determine the slip rate for a given stress and strength value. In addition, we used an aging law16,34, which can be considered an evolution law for changes in strength that varies depending on prior slip history. In our simulation of earthquake generation cycles, we used the quasi-dynamic (QD) model32 with a smaller radiation damping term instead of a fully dynamic (FD) model. Usually, the QD model cannot simulate fast slip during seismic events; however, in models with standard rate-and-state friction and relatively uniform fault properties, the FD and QD simulations produce qualitatively similar fault behaviors, with crack-like ruptures and similar earthquake patterns35.

During our simulations, we only changed the numerical values and the spatial distributions of the frictional parameters (A, B, and L), focusing on the slip in the shallow, southern area. Nakata et al.31 approximated the spatial distributions of frictional parameters for seismic sources with circular patches for M7 earthquakes, which repeatedly occurred at intervals of several decades in the past. In the present study, we removed these circular patches for simplicity and used the model without these unstable patches (Model B1, Supplementary Fig. 3, Supplementary Table 2), instead of Nakata et al.31, for comparison. The calculation time was then reduced so that we could calculate multiple cycles of the M9 earthquakes. In a method similar to that of Nakata et al.31, the M9 coseismic rupture area located in the middle segment of the Japan Trench was approximated by a rectangle (an area 150 km long along the strike and shallower than 22 km in depth). At the boundary between the southern and northern sides of the M9 source area, the frictional characteristics changed gradually, and differences between the deep and shallow areas were not considered.

For Model B3 (Supplementary Fig. 4, Supplementary Table 3), the values and spatial distribution of the frictional parameters (A, B, and L) were determined with reference to observations of the residual gravity anomaly by Bassett et al.14. The M9 coseismic rupture area was approximated by a rectangle (170 km along the strike and shallower than 24 km depth). Here, we divided the surrounding area, which we assumed to have a uniform frictional condition as in Model B1, into three segments: northern, southern, and deep segments (Supplementary Fig. 4). Difference of frictional parameters at the boundaries between the northern, southern, and middle (the M9 coseismic) segments were then assumed to be larger than those of Model B1 (Supplementary Fig. 3c–d), and their locations were based on an observed residual gravity anomaly < 0. Spatial distributions of the coseismic slip, postseismic slip, seismic velocity, and residual gravity anomaly showed that uniform and smooth assumptions of frictional boundaries may not be appropriate. Boundary settings along the trench between the middle and southern segments were the main differences between Model B1 and Model B3. In addition, based on our hypothesis, the frictional parameter controlling slip weakening (L) at the southern segment is assumed to be larger than that of the middle and deep segments. Through iterative processes of trial and error, fine-tuning of the frictional parameters was undertaken between Models B1 and B3.

To assess the effect of the channel layer, we conducted calculations using four models (Supplementary Table 1). Model B1 was almost the same as in our previous study. This model did not have a channel layer. Models B2, B3, and B4 were the modified models in this study. Model B2 had a small L value (the same as Model B1), Model B3 had a medium L value, and Model B4 had a large L value in the southern segment. Model B2 had a channel layer in which contrast between the middle and southern segments was weak. Models B3 and B4 represented models with a high-contrast channel layer. The L values in the southern segment were 2.1, 4.6, and 6.4 times larger than in the middle segment for Models B2, B3, and B4, respectively.