Nakata et al.
Progress in Earth and Planetary Science(2023) 10:34
https://doi.org/10.1186/s40645-023-00566-y

RESEARCH ARTICLE

Progress in Earth andPlanetary Science

Open Access

Recurrence intervals for M > 7
Miyagi‑ken‑Oki earthquakes during an M ~ 9earthquake cycle
Ryoko Nakata1* , Naoki Uchida2,3, Takane Hori4 and Ryota Hino2

Abstract
The 2011 Tohoku-Oki great earthquake increased the difficulty of evaluating the long-term probability of seismic
activity along the Japan Trench because of the unknown impact of the unprecedentedly large slip. In this study, thelong-term activity of “Miyagi-ken-Oki earthquakes”, an M > 7 earthquake sequence off Miyagi Prefecture, located at
the edge of the source area of the Tohoku-Oki earthquake was simulated. We conducted numerical simulations ofearthquake generation cycles based on the rate- and state-dependent friction law representing the stress accumulation and release process on the plate interface. We also validated the results based on analyses of observationaldata, including time intervals and afterslip distributions for the previous Miyagi-ken-Oki earthquakes. The simulation
results were then compared with repeating-earthquake-derived interplate slip observations over 30 years. The results
showed that the spatial and temporal characteristics of aseismic slips on the plate interface near the M > 7 Miyagiken-Oki source were qualitatively reproduced. The time interval between the M ~ 9 earthquake and the first M > 7
earthquake is shorter than the average recurrence interval of M > 7 earthquakes during the latter stage of each M ~ 9earthquake cycle. In contrast, the interval between the first and the second M > 7 earthquakes is the longest in eachM ~ 9 earthquake cycle. The time intervals of the M > 7 earthquakes fluctuated in the early stage compared to thosein the latter stage of the M ~ 9 earthquake cycle. These characteristics were associated with the extent of the locked
and afterslip areas in and around the source. Hence, monitoring the spatio-temporal distribution of the slip rate in
and around the source area during the preparation process of earthquakes occurring in the downdip area off Miyagi
Prefecture is very important to assess whether the next M > 7 earthquake is approaching. Furthermore, earthquake
cycle simulations combined with earthquake and slow slip monitoring can provide meaningful insights for long-term
assessments of great interplate earthquakes.
Keywords Earthquake cycle, Recurrence interval, Slip deficit, Tohoku-Oki earthquake, Miyagi-ken-Oki earthquake,
Repeating earthquake, Aseismic slip, Rate- and state-dependent friction law*Correspondence:
Ryoko Nakata
rnakata@eps.s.u-tokyo.ac.jp
Full list of author information is available at the end of the article© The Author(s) 2023. Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which
permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to theoriginal author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images orother third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line
to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory
regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of thislicence, visit http://​creat​iveco​mmons.​org/​licen​ses/​by/4.​0/.

Nakata et al. Progress in Earth and Planetary Science(2023) 10:34

Page 2 of 16

Fig. 1 Distribution of the source areas of recent large interplate earthquakes off Miyagi prefecture and the 2011 Tohoku-Oki earthquake. Blackdashed contours: coseismic slip distribution (10 and 50 m) of the 2011 Tohoku-Oki earthquake (Iinuma et al. 2012). Gray lines: distribution of theestimated postseismic slip (> 0.4 m) for 0.6 years following the 2011 Tohoku-Oki earthquake (Iinuma et al. 2016). Black and gray stars indicate thehypocenters of the mainshock (Mj = 9.0) and foreshock (Mj = 7.3) of the 2011 event, respectively, determined by JMA. The magenta and orange
areas indicate the source areas (slip > 0.6 m) of the 1978 and 2005 Miyagi-ken-Oki earthquakes, respectively (Yamanaka and Kikuchi 2004; Yamanaka
2005). The red and green stars indicate the hypocenters of the March and May 2021 earthquakes determined by JMA. The gray region on landsurrounded by the solid line indicates the location of Miyagi Prefecture. The brown line indicates the plate boundary (Bird 2003). The rectangle in
the inset with the plate boundary (Bird 2003) represents the location of the study area. PA Pacific plate, PS Philippine Sea Plate, AM Amur plate, OK
Okhotsk plate

1 Introduction
In the middle segment on the west side of the Japan
Trench, M ~ 9 and M > 7 interplate earthquakes haveoccurred on the updip and downdip parts of the subducting Pacific Plate interface, respectively (Fig. 1). The
2011 Tohoku-Oki earthquake at Mj9.0 (Mj refers to themagnitude determined by the Japan Meteorological
Agency (JMA)) occurred on the updip side of the segment. Geological data, tsunami sediment, and historical documents indicate that the recurrence interval of
great tsunamigenic earthquakes along the Japan Trenchis < 1000 years (Uchida and Burgmann 2021). Whereas,

Nakata et al. Progress in Earth and Planetary Science(2023) 10:34

the downdip edge of the coseismic slip area in the 2011
Tohoku-Oki earthquake was located off the east coastof Miyagi Prefecture, where large (Mj7.1–7.4) earthquakes have occurred several times in the past, particularly in 1897 (M7.4) (Utsu 1979), the 1930s (1933 (Mj7.1),
1936 (Mj7.4), and 1937 (Mj7.1)), 1978 (Mj7.4), and 2005
(Mj7.2) (Headquarters for Earthquake Research Promotion 2011). The 1930s earthquakes (1933, 1936, and 1937)
are treated as a single event series based on the source
regions. The rupture areas of these earthquakes, which
are referred to as “Miyagi-ken-Oki earthquakes,” partially overlapped, with maximum coseismic slips < 3 m(Yamanaka and Kikuchi 2004; Yamanaka 2005).
The probability of an earthquake occurring within
30 years is the most commonly used indicator for longterm evaluations provided by the Headquarters for
Earthquake Research Promotion (HERP), a Japanese
government organization responsible for earthquake
research. The probability is calculated basically assuming a BPT (Browning Passage Time) distribution model(Matthews et al. 2002) using the occurrence time of
the most recent earthquake and the mean recurrence
interval of the earthquakes in the region, while Poissonprocess is assumed if the occurrence time of the mostrecent event is unknown. For the evaluation of the Miyagi-ken-Oki earthquake by HERP, the 2011 Tohoku-Oki
earthquake was treated as the most recent earthquake.
This is because the Tohoku-Oki earthquake ruptured
a wide region including the source regions of the previous Miyagi-ken-Oki earthquakes. The slip amount of
the Tohoku-Oki earthquake was estimated to be ~ 10 min the Miyagi-ken-Oki source area (Iinuma et al. 2012)(Fig. 1), which is larger than those in previous Miyagiken-Oki events (Nakata et al. 2016). The mean recurrence
interval was calculated to be 38.0 years using 1897, 1936,
1978, and 2011 events. The 2005 event was not used in
the HERP’s evaluation because the event ruptured only apart of the previous Miyagi-ken-Oki earthquake sourcearea (Yaginuma et al. 2007). Applying the average intervalof 38.0 years and elapsed time of 10.8 years (approximate
as of January 1, 2022) to the BPT distribution model, theprobability of an earthquake within the next 30 years is
estimated to be ~ 70–80% (HERP 2022).
Approximately 10 years after the Tohoku-Oki earthquake, Mj6.9 and 6.8 interplate earthquakes occurred in
the region off Miyagi Prefecture on March 20 and May 1,
2021, respectively (Fig. 1). The coseismic slip area of thetwo earthquakes comprised only the western part of thesource region for the 1978 Miyagi-ken-Oki earthquake(Mj7.4) (HERP 2021a, 2021b; Yoshida et al. 2022), which
was one of the largest earthquakes in the Miyagi-ken-Oki
segment. Although the 2021 earthquakes were smaller
in magnitude than the ~ M7.4 earthquake predicted toPage 3 of 16

occur off Miyagi Prefecture (~ 70–80%) within the next
30 years by HERP, they were the largest interplate earthquakes to rupture within the 1978 source area since the2011 Tohoku-Oki earthquake.
Nakata et al. (2016) investigated the seismicity after agreat (M ~ 9) earthquake using three-dimensional numerical simulations. They suggested that the M > 7 Miyagiken-Oki (MYG) earthquake after the great earthquake
would occur earlier than the average recurrence interval
during the later stage of the interseismic period of thegreat earthquake. This is due to the locking conditions
around the MYG patch and high-stress rate caused by
the large afterslip of the M ~ 9 event at the downdip sideof the M ~ 9 rupture area, promoting earthquake occurrence (Nakata et al. 2016). Although Nakata et al. (2016)
have modeled five M ~ 7 earthquakes off Miyagi (MYG),
off Fukushima (FKS), off Ibaraki (IBK), and at the shallow
regions (SHL1 and SHL2), comparisons with observations were made qualitatively for recurrence interval andmagnitude of each earthquake that occurred prior to the2011 Tohoku-Oki earthquake.
Although Nakata et al. (2016) only focused on MYG
earthquakes that occurred in the 200 years prior to an
M ~ 9 earthquake and one MYG earthquake after theM ~ 9 earthquake, to reliably propose the long-term seismic activity in the MYG region, it is necessary to investigate the entire M ~ 9 earthquake cycle after knowing how
well the simulation results reasonably account for recent
observations. However, a long time was required to calculate multiple M ~ 9 cycles using a model with manyM ~ 7 patches. Furthermore, the propagation of afterslip
due to earthquakes in other areas complicates the study
and description of the relationship between the MYG
earthquakes and the M ~ 9 source, despite other earthquakes not being expected to significantly affect the timeinterval of MYG earthquakes in the earliest stage of theM9 cycle (Nakata et al. 2016). We then prepared new
simulation results to discuss the recurrence interval of
MYG earthquakes through the M9 earthquake cycle.
Accordingly, in this study, we examined detailed spatial and temporal characteristics of simulated slips on theplate interface off Miyagi Prefecture common to multiple
great earthquake cycles obtained from one model along
the Japan Trench. We then validated our findings by evaluating how well the simulation results reproduced observations since 2011. Using a similar approach as in ourprevious studies along the Japan Trench (Nakata et al.
2016, 2021), we numerically simulated earthquake generation cycles based on stress accumulation and releaseprocesses on the plate interface. The work in Nakata et al.
(2021) is a revision of that in Nakata et al. (2016) withrespect to the numerical and spatial settings of friction
parameters along the trench to improve reproducibility

Nakata et al. Progress in Earth and Planetary Science(2023) 10:34

for the observed M ~ 9 source area and long afterslip
duration. Since Nakata et al. (2021) does not account for
the frictional heterogeneity that represents M ~ 7 earthquakes, including those of the Miyagi-ken-Oki earthquake, we performed a new calculation combining thefollowing: the M ~ 9 and Miyagi-ken-Oki earthquakes.
Hence, the results of this study will improve our understanding of the basic physical processes associated withearthquake occurrence and evolution in the Miyagi-kenOki region throughout M ~ 9 earthquake cycles.

2 Methods
The equations, initial conditions, basic parameters, anddiscretized plate geometry were the same as those usedin previous studies (Nakata et al. 2016, 2021), as briefly
introduced below. Details are provided in Additional
file 1.
2.1 Numerical simulation of earthquake generation cycles

We discretized the subducting Pacific Plate interface
into small planar triangular elements (Hyodo et al. 2016;
Additional file 1: 1.2–1.3) and treated three triangularelements as a rectangle-like subfault (Additional file 1:
Fig. S1). Seismic and aseismic events release the slip deficit accumulated during interseismic periods. The time
derivative of shear stress on a rectangle-like subfault
caused by slip rate changes of surrounding subfaults wasevaluated by the combined effect of angular dislocations
by three triangular elements within a rectangle-like subfault (Comninou and Dundurs 1975) in a homogenous
elastic half-space (Hyodo et al. 2016). The quasi-dynamic
model (Rice 1993) was used to treat the long cycle. We
employed a rate- and state-dependent frictional law (Dieterich 1979) as an approximate mathematical model for
friction behavior at the plate interface and a fault constitutive law (Nakatani 2001) to determine the interplate
slip rate for a stress and strength value. We also employedan aging law (Dieterich 1979; Ruina 1983). The frictionalparameters were assumed to be time-invariant in thesimulations (Nakata et al. 2016, 2021).
2.2 Friction parameters

We assumed friction parameters A, B, and L for each subfault (Fig. 2) based on those described by Nakata et al.
(2016, 2021) to satisfy the following conditions for nondimensional numbers Rb and Ru (Additional file 1: 1.4),
which controlled the characteristics of the simulated slip
(Barbot 2020). The friction parameter A (= aσ) controls
the slip increase rate at which the level of stress increases
and reaches the strength; friction parameters B (= bσ)
and L control the strength recovery and slip weakening (Additional file 1: 1.1). When Ru > 1, the fault slip isunstable (Barbot 2020); however, when Ru ~ 1 and < 1, thePage 4 of 16

fault slip is quasi-stable, and stable, respectively (Barbot2020). When the Rb value is small and positive, conditions approach a neutral velocity; whereas when the Rb
value is large and approaches 1, conditions correspond toa notable weakening unless a is unusually large (Barbot2020).
The friction parameters in the aseismic backgroundarea (northern, southern, and deep segments) were thesame as those described for Model B3 in Nakata et al.
(2021). The M ~ 9 source area assumed in the middle segment was expressed by a rectangle-like region (Nakata
et al. 2021). We expressed the source area of the historical M > 7 Miyagi-ken-Oki (MYG) earthquakes using acircular patch when projected onto a horizontal plane(Nakata et al. 2016). Although the spatial distribution was
the same, we slightly changed the friction parameter (A,B, and L) values for these seismic sources from those in
our previous studies (Nakata et al. 2016, 2021). In addition to A–B and L, the above-mentioned nondimensionalparameters are shown in Table 1. In these seismic areas,
Ru > 1 with a low Rb.

3 Results
We performed extensive numerical simulations forapproximately 3300 years. Figure 3 shows the development over time of the slip deficit obtained using themodel shown in Table 1 and Fig. 2. Here, we defined aslip deficit as Vpl × time − (slip at the point) + α, where Vpl
is the plate convergence rate (Additional file 1: Table S1)and α is a value-added for visibility and does not changewith time. In Fig. 3, the slip deficit at the cross in theM ~ 9 slip area (shallow middle segment in Fig. 2a) was
set to zero by adjusting α when the first M ~ 9 earthquakeof Cycle 3 occurred (T = 1919.3 years). At that point, fivelarge > 50 m changes (Fig. 3a, black line) caused by M ~ 9(M = 8.94–8.95) earthquakes (such as Fig. 4a in the slip
distribution) were observed, with recurrence intervalsof 633–637 years. These periods between two successive
M ~ 9 earthquakes are hereafter referred to as “M9 earthquake cycles.”The slip deficit at the cross in the patch labeled “MYG”
in Fig. 2a was set to zero when the M > 7 event immediately following the first M ~ 9 earthquake of Cycle 3
occurred (T = 1968.7 years in Fig. 3a). The time development of the slip deficit at the cross in the MYG patchat the downdip region of the middle segment (Fig. 3,
magenta line) revealed the occurrence of numerous
earthquakes with a slip of approximately 5 m. Over a
period of 3300 years, earthquakes with magnitudes of
7.34–7.45 (such as those in Figs. 4b and 4c presenting
the slip distributions) repeatedly occurred in the MYG
patch at time intervals of 60–123 years (Figs. 3 and 5).
As the source area was approximated by a single circle

Nakata et al. Progress in Earth and Planetary Science(2023) 10:34

Page 5 of 16

Fig. 2 Distribution of frictional parameters used in this study. a Spatial distribution on the plate interface of (A–B) (MPa). Contours: depth (km) to
the upper surface of the descending plate (Baba et al. 2006). Gray solid curve: downdip limit of the interplate earthquake distribution (Igarashi et al.
2001). Black crosses: locations used in Fig. 3. Barbed line: trench axis. b Spatial distribution of the characteristic slip distance (L). c–d Profiles of A − B
and L along the strike at a depth of 10.15 km. MYG, Miyagi-ken-Oki

and homogeneous frictional conditions were assumed
within the circle, minimal fluctuation in the magnitudeof each earthquake during the M9 earthquake cycle was
expected. However, the maximum slip and magnitudeof the MYG earthquakes tended to be slightly smaller
immediately after the M ~ 9 earthquake (Figs. 4c, 5e, andn) than those of the other earthquakes (Figs. 4b, 5a–c,
and f–l). In addition, the coseismic slip areas of the MYG
earthquakes were slightly narrower in dip direction due
to a lack of slip on the updip side immediately after theM ~ 9 earthquake (Figs. 4c and 6, red lines) when compared to those of the other MYG earthquakes (Figs. 4band 6, blue lines), regardless of the magnitude of each
MYG earthquake.
Notably, immediately after the M ~ 9 earthquake, theresulting aseismic slip distribution after the M > 7 MYG
earthquake also differed from that of the other earthquakes. During the M9 earthquake cycle, aseismic slip
was distributed around the patch and sometimes widely
propagated to the updip (Fig. 6, light-blue lines). In

contrast, the aseismic slip did not typically propagate to
the updip side immediately after the M ~ 9 earthquake(Fig. 6, pink lines).
The recurrence intervals between MYG earthquakes
showed lengthening and shortening periods during theM9 earthquake cycles. The patterns of lengthening or
shortening were different in Cycles 1 to 4 in the latter halfof the M9 earthquake cycle. Among them, the common
characteristics were as follows. At the latest phase of each
M9 earthquake cycle, the average recurrence intervalof the four M > 7 seismic events in the MYG patch was
69–70 years. In contrast, the variability of the recurrence
interval of the MYG earthquake was large during an initial phase of the M9 earthquake cycle. The time interval
between each M ~ 9 earthquake and subsequent M > 7
MYG earthquake was 49–57 years, which were shorter
than the average recurrence intervals of the MYG earthquakes during each M9 earthquake cycle. Furthermore,
the next MYG event had the longest interval from theprevious one (> 100 years).

Nakata et al. Progress in Earth and Planetary Science(2023) 10:34

Page 6 of 16

Table 1 Parameters at the northern, middle, southern, and deep segments and the MYG patch
Segments
Northern

MiddleSouthern

Parameters

Values

Sources/Segments

A–B (MPa)− 0.10MYGb–a

0.0020Parameters

ValuesA–B (MPa)− 0.320b–a0.0064

L (m)

0.60

L (m)

0.026

Rb = (b–a)/b

0.118

Rb = (b–a)/b

0.299

Ru

0.271

Ru

14.2

Length (km)*2

65 (*1)Radius (km)23

A–B (MPa)− 0.19b–a0.0038

Deep (≤ 62 km *3)

A–B (MPa)− 0.12b–a0.0024

L (m)

0.16

L (m)

0.15

Rb = (b–a)/b

0.202

Rb = (b–a)/b

0.138

Ru

5.047

Ru

9.6

Length (km)*2

170

Length (km)*2

480

A–B (MPa)− 0.10A–B (MPa)

0.10

b–a

0.0020b–a− 0.0020

L (m)

13.0

Length (km)*2

480

L (m)

0.65

Rb = (b–a)/b

0.118

Ru

0.865

Length (km)*2

225 (*1)

Deep (> 62 km *3)*

1: Transition between the middle segment and northern/southern segment is 10 km*

2: Length along strike direction*

3: Changes depending on the down-dip limit of the interplate earthquake distribution (Igarashi et al. 2001)

Fig. 3 Temporal distribution of slip deficits obtained using the friction model shown in Fig. 2. a Temporal distribution of simulated slip deficits overa 3300-year period at a point within the M ~ 9 source area (black) and Miyagi-ken-Oki (MYG) source area (magenta) at the locations of the crosses
shown in Fig. 2. Cycles 1–4 (light gray lines) represent a period from one M ~ 9 earthquake to the next M ~ 9 earthquake. The slip deficit in the M ~ 9
area was set to zero when the first M ~ 9 earthquake of Cycle 3 occurred (T = 1919.3 years), and the vertical axis is on the left. The slip deficit at the
MYG area was set to zero when the MYG earthquake occurred immediately after the first M ~ 9 earthquake of Cycle 3 (T = 1968.7 years), and thevertical axis is on the right. b Temporal distribution of the simulated slip deficit at the cross in the MYG patch shown in Fig. 2a for a 1000-year period
from 200 years before the occurrence time of each M ~ 9 earthquake in Cycles 1–4. The time when the slip deficit is zero differs from that in (a), andthe slip deficit is a relative value for ease of viewing. The period combines the dark and light gray lines in (a). Horizontal axis (T′): 0 years at the timeof the M ~ 9 earthquake, which is the start of each cycle. c Enlarged diagram of (b) for T′ = − 15 to 15 years

Nakata et al. Progress in Earth and Planetary Science(2023) 10:34

Page 7 of 16

Fig. 4 Coseismic slip distribution obtained using the friction model shown in Fig. 2. a Simulated coseismic slip distribution of the M = 8.95
earthquake that occurred at T′ = 0 years in Cycle 3. Magenta and white crosses: locations used in Fig. 3. The large square surrounding the magenta
cross indicates the area described in (b) and (c). b Simulated coseismic slip distribution for one of the events in the Miyagi-ken-Oki (MYG) patch in
the last stage of the M ~ 9 earthquake cycle, the M = 7.34 earthquake at T′ = − 1.7 years in Cycle 3. Magenta cross: location used in Fig. 3. c Simulatedcoseismic slip distribution for one of the first events after an M ~ 9 earthquake, the M = 7.35 earthquake at T′ = 49.4 years in Cycle 3

At the end of each M9 earthquake cycle, the interval
between the last MYG earthquake and M ~ 9 earthquake
was 0.8–2.1 years (Fig. 3c), which was much shorter
than the average recurrence interval of the MYG events;on the other hand, the interval and magnitude of earthquakes occurring in the MYG patch at the final stage ofeach M9 earthquake cycle varied between cycles. These
characteristics were common to all cycles over the evaluated 3300-year period (Fig. 3b). Although a small slow
slip event occurred near the MYG patch towards the endof Cycle 2 (~ 550 years in Fig. 3b), this phenomenon did
not affect the overall recurrence interval and aseismic
slip distribution trends.

4 Discussion
During an M ~ 9 earthquake cycle, the stress and strain
at the plate boundary vary spatiotemporally, and thesecontinuous variations are expected to cause changes inthe time interval of MYG earthquakes. In other words,
changes in recurrence intervals of MYG earthquakes
may reflect such changes in physical conditions. In thissection, we divided an M ~ 9 earthquake cycle into threestages and extracted characteristics common to all M ~ 9earthquake cycles: from an M ~ 9 earthquake to thesubsequent MYG earthquake (Stage I: earliest stage),from the first MYG to the next MYG earthquake (Stage
II: early stage), and from the second MYG to the nextM ~ 9 earthquake (Stage III: later stage). We first compared the simulation results with observational studies
to assess the validity of the simulations. We then identified factors influencing the recurrence intervals of M > 7

earthquakes at each stage. To this end, the results were
examined using multiple patch models with other M ~ 7
earthquakes at the updip and southern regions and thecorresponding stage to the post-2021 earthquakes off
Miyagi among the three stages. Finally, we presented
several important points for the long-term evaluation of
M > 7 earthquakes off Miyagi Prefecture.
4.1 Comparison between observations and simulations

Based on more than 110 years of observations since 1897,
the variability in the recurrence intervals of the Miyagiken-Oki earthquakes has been relatively small. In thesimulation, at the later stage of the M9 earthquake cycle(Stage III), the average of the recurrence intervals for the
M > 7 events in the MYG patch was relatively constant. At
the end of each M9 earthquake cycle, however, the interval between the last MYG event and the next M ~ 9 earthquake in the simulations was significantly shorter thanthe observed interval between the 1978 and 2011 events.
The interval between the 2005 Miyagi-ken-Oki earthquake and the 2011 Tohoku-Oki earthquake, on theother hand, was much shorter than the observed average recurrence interval of Miyagi-ken-Oki earthquakes.
This short interval may be related to the extremely shorttime between the last MYG earthquake and the nextM ~ 9 event in each M9 earthquake cycle obtained in thesimulations. Afterslip of the last MYG earthquake might
have facilitated the occurrence of the subsequent M ~ 9earthquake, as suggested for the 2005 Miyagi-ken-Oki
earthquake (Suito et al. 2011; Ozawa et al. 2012). If so,
it may have significant implications for predicting theNakata et al. Progress in Earth and Planetary Science(2023) 10:34

Page 8 of 16

Fig. 5 Snapshots of the simulated slip velocity distribution during Cycle 3 shown in Fig. 3. a–c M > 7 earthquake in the Miyagi-ken-Oki (MYG) patch.
Slip velocity is normalized to the plate convergence rate, Vpl. The time interval between adjacent panels is shown below the arrow between thepanels. d and m M ~ 9 earthquakes that occurred within a time interval of 633 years. e–l, n M > 7 earthquakes that repeatedly occurred in the MYG
patch. Blue and red areas indicate the locked parts and unstable slipping parts of the fault, respectively. Yellow/green and white indicate slow slip
and plate convergence rates, respectively

occurrence of great interplate earthquakes, but we will
not go into this issue further in this study.
Differences in slip distribution described by the simulations in this study were consistent with the observations from past earthquakes. Widely distributed
aseismic slips during Stage III were consistent with theobserved aseismic slips of the 2005 and 1978 earthquakes. The aseismic slip extending to the updip side
was evidenced by the onshore geodetic observation(Miura et al. 2007) and the increased activity of repeating earthquakes (Uchida et al. 2006, 2007) after the2005 earthquake. Although geodetic observations only
detected the deep aseismic slip after the 1978 earthquake, Uchida et al. (2006) and Hino et al. (2007)
pointed out remarkable similarity between the aftershock distributions of the 1978 and 2005 earthquakes
on the updip side of the mainshock ruptures. The 2005
aftershocks included repeating earthquakes, suggesting

Nakata et al. Progress in Earth and Planetary Science(2023) 10:34

Page 9 of 16

Fig. 6 Coseismic and aseismic slip distributions of M > 7 earthquakes around the Miyagi-ken-Oki (MYG) patch. Coseismic slip is extracted whenslip velocity (V) > 1 cm/year. Aseismic slip is the slip remaining after removing coseismic slip from cumulative slip over the three years from therupture initiation time of each M > 7 earthquake. a Map view of simulated coseismic slip (> 1 m) and aseismic slip (> 1 m) distributions of the MYG
earthquakes shown in Fig. 5. The coseismic slip of the MYG earthquakes before the M ~ 9 earthquake (10 events; Fig. 5a–c, f–l) is indicated in blue,
and light-blue lines show the aseismic slip following these earthquakes. Among the light-blue lines, the dotted lines indicate events that occurredjust before M ~ 9 earthquakes with only one year of cumulative slip (Fig. 5c and l). The reason for the one-year duration is owing to the shortest timeinterval of < 3 years between the MYG earthquake and M ~ 9 earthquake. The coseismic slip of the MYG earthquakes after the M ~ 9 earthquake (2events; Fig. 5e and n) is indicated in red, and pink lines show the aseismic slip following these MYG earthquakes. Black crosses represent the points
shown in Fig. 7. b and c Slip profiles of coseismic slip and aseismic slip distributions of the MYG earthquakes shown in Fig. 5 along gray lines (C–C′,
across the MYG patch center; F–F′, along 35 km depth) in the map panelthat the shallow 1978 aftershocks were due to aseismic
slip as observed in the 2005 activity. The slightly narrower coseismic slip distributions simulated in Stage I
could be considered similar to those of the 2021 earthquakes (Yoshida et al. 2022), which ruptured exclusively
on the downdip side of the rupture areas of earthquakes

in 1978 (Yamanaka and Kikuchi 2004) and the 1930s(Umino et al. 2007).
To perform comparisons with observations of small
to medium earthquakes, which could not be addressedin this simulation, cumulative slips are shown (Fig. 7).
The development over time of the averaged cumulativeNakata et al. Progress in Earth and Planetary Science(2023) 10:34

Page 10 of 16

Fig. 7 Cumulative slip in the middle segment before and after an M ~ 9 earthquake. a Cumulative slip from 1993–2021 estimated from observed
repeating earthquake groups within each rectangle. Vertical line: time of 2005 off Miyagi and 2011 Tohoku-Oki earthquakes. b Location of eachrectangle. Crosses: locations shown in (c). c Simulated cumulative slip time series over 30 years (from T′ = − 19 years to T′ = 11 years) in Cycle 3.
Points at crosses in regions a1, a2, a3, b1, b2, and c1 shown in (b) were used. The behavior within each rectangle was represented by the points
marked with crosses. Solid lines: from T′ = − 19 to 0 years. Dotted lines: from T′ = 0.2 to 11 years by subtracting the slip due to the M ~ 9 earthquake.
The horizontal axis (T′): 0 years at the time of the M ~ 9 earthquake. Dashed gray lines: slope at Vpl. Vertical lines: time of M > 7 and M ~ 9 earthquakes.
d Enlarged view of (c) at point b1. MYG, Miyagi-ken-Oki

interplate slip estimated from repeating earthquakes(Uchida and Matsuzawa 2013) exhibited a heterogeneous pattern. As of 2021, the curves could be categorized into three groups: regions showing large slip after
the Tohoku-Oki earthquake (group a; a1, a2, and a3 inFig. 7a–b), small slip after the earthquake (group b; b1
and b2), and no postseismic slip (group c; c1). Moreover, following the Tohoku-Oki earthquake, many repeating earthquakes emerged on the downdip side that wereaseismic before the Tohoku-Oki earthquake while interplate earthquakes were relatively absent from inside
the Tohoku-Oki rupture area (Yoshida et al. 2022). The
former corresponded to group a, while the latter corresponded to group c. Similar results were observed for thesimulation after an M ~ 9 earthquake. That is, immediately after an M ~ 9 earthquake (Stage I), a large afterslip
on the deeper side off Miyagi Prefecture (a1, a2, and a3
in Fig. 7c) was reproduced. Furthermore, a small aseismic

slip occurred between the area off Miyagi Prefecture andthe M ~ 9 area (b1 and b2 in Fig. 7c and d). A point on
the updip side (c1 in Fig. 7c) became nearly locked, whichis reasonable given that c1 is located within the unstablearea ruptured in the M ~ 9 earthquake (Additional file 2:
Fig. S2.1). Although the amount of slip varied, after an
M ~ 9 earthquake, the slip on the updip side (b1, b2, andc1) was sufficiently small compared to that on the downdip side (a1, a2, and a3). These results were consistent
with prior observations (Fig. 7a).
Observations of repeating earthquakes showed a small
increase in slip in area b1 after the 2005 earthquake(M7.2), with no change detected in area a2 (Fig. 7a). In
contrast, geodetic data revealed a small afterslip aroundthe a2 area (Miura et al. 2007). In the simulation results(Figs. 6a and 7c), the slips clearly increased at points
a2, b1, and b2 after the MYG earthquake (M7.34) that
occurred ~ 1.7 years before the M ~ 9 earthquake.

Nakata et al. Progress in Earth and Planetary Science(2023) 10:34

Page 11 of 16

Fig. 8 a Time development of simulated slip velocity during Cycle 3. Horizontal axis: depth (km) through the approximate center of the
Miyagi-ken-Oki (MYG) patch (along the C–C′ line). Vertical axis: 0 years at the time of the first M ~ 9 earthquake during Cycle 3 (similar to Fig. 3). Red
horizontal lines: timing of M > 7 or M ~ 9 earthquakes shown in Fig. 5. Black dashed line: center of the MYG patch; solid vertical lines: a radius of
23 km from the center of MYG patch. b Enlarged diagram of (a) for T′ = − 15 to 15 years

The amount of simulated aseismic slips was larger
than that in the observations. However, the estimatedmoment during 2.72 years for the postseismic slip of the2005 earthquake (12.19 × ­1019 Nm in Suito et al. 2011)
was comparable with the three-year aseismic slip estimated at the later stage of the M9 cycle in the simulation(8.7–11.6 × ­1019 Nm). Therefore, the simulation reproduced the overall seismicity features in the Miyagi-kenOki region, although the simulation and observations
sometimes differed with respect to factors such as thetime intervals between each earthquake and aseismic slip
amounts.
4.2 Factors influencing recurrence intervals of M > 7
Miyagi‑ken‑Oki earthquakes

In Sect. 4.1, we showed that the simulation results in thisstudy were sufficiently reasonable to describe the difference in the spatiotemporal behavior of slips (both seismic
and aseismic) before and after an M ~ 9 earthquake whencompared to prior observations. Although the descriptions might not be sufficient, they provide a necessary
comparison to validate our simulation results. Therefore,

we investigated the physical process responsible for thedifferences in the recurrence intervals among the threestages.
The trend in the time evolution of the slip deficit in
the center of the MYG patch (Fig. 3, magenta line) was
variable throughout an M9 earthquake cycle. More specifically, although the change in the slip deficit was large
in Stages I and II, the change in the slip deficit was relatively low and nearly constant in Stage III as the nextM ~ 9 earthquake approached. We then compared thespatial and temporal changes in slip velocity and stresswithin and around the MYG patch, along dip and strikedirections in each stage (Figs. 8 and 9, Additional file 2:
Figs. S2.2–S2.5).
During Stage III (~ T’ = -200–0 and 170–630 years in
Fig. 8), immediately after each MYG earthquake, the slip
velocity within the patch became low (locked condition)and the afterslip was distributed both in the updip anddowndip direction around the patch. The locked area
within the patch (Fig. 8 and Additional file 2: Figs. S2.3–
S2.4) gradually narrowed from the outer edge of thepatch toward the center as the next M > 7 earthquake

Nakata et al. Progress in Earth and Planetary Science(2023) 10:34

Page 12 of 16

Fig. 9 a Time development of simulated stress during Cycle 3. Horizontal axis: depth (km) through the approximate center of the Miyagi-ken-Oki(MYG) patch (along the C–C′ line). Vertical axis: 0 years at the time of the first M ~ 9 earthquake during Cycle 3 (similar to Fig. 3). Red horizontal lines:
timing of M > 7 or M ~ 9 earthquakes shown in Fig. 5. Black dashed line: center of the MYG patch; solid vertical lines: a radius of 23 km from thecenter of MYG patch. b Enlarged diagram of (a) for T′ = − 15 to 15 years

approached, whereas, the center of the patch remained
locked. Eventually, the inside of the patch became unstable, slow slip occurred around the edges within the patch,and fast slip (i.e., an earthquake) occurred later.
The difference in the slip distribution and recurrence
interval between Stage I (~ T’ = 0–50 and 630–690 years
in Fig. 8) and the other stages might be attributed to thewide distribution of the locking area due to the M ~ 9earthquake. After an M ~ 9 earthquake, due to the stronglocking of the updip side of the MYG patch (Figs. 5e,5n, and 8), the afterslip following the first MYG earthquake did not propagate substantially to the updip side
(Figs. 6 and 8). The afterslip following an M ~ 9 earthquake showed a higher slip velocity (Nakata et al. 2016)
and longer duration only on the downdip side than in theother stages. In addition, the afterslip penetrated inside
the patch, and the locked area within the patch was narrower than in the other stages. The stress level in thestrongly locked patch was higher in Stage I than in theother stages (Fig. 9). Finally, the first MYG earthquake
occurred earlier due to the high stress accumulation rate
(Nakata et al. 2016).

During Stage II (~ T’ = 50–170 years in Fig. 8), the locking on the downdip side of the patch was similar to thatin Stage III. Given that the afterslip area was especiallynarrow and weak on the updip side of the MYG patch,the locked area on the updip edge (approximately 27 km
depth) remained wider than in Stage III. On the updip
side of the locked area within the patch, it subsequently
narrowed more slowly than that in Stage III. Stress in thepatch remained low during the early period of Stage II(Fig. 9). Therefore, the next M > 7 earthquake will occur
following a longer time interval than that in Stage III, as
the interior of the patch will remain relatively stable for alonger period.
In numerical simulations using a two-dimensionalrheological model at the Japan Trench (Barbot 2020),
M ~ 7 earthquakes occurred with a recurrence intervalthat increased during an M ~ 9 earthquake cycle. Simply
stated, the recurrence interval was longer in the later partof the M9 earthquake cycle than in the beginning. This
simple situation can be considered similar to this study,
where Stage I has a shorter recurrence interval than Stage
III. However, the behavior at the area between the M ~ 9

Nakata et al. Progress in Earth and Planetary Science(2023) 10:34

large slip area and M ~ 7 earthquake during Stages I and
II differed from the current results. More specifically,in the current study simulation, this area was locked, oran exceptionally small slip occurred between the areas(Figs. 7 and 8). In contrast, in the simulations conductedby Barbot (2020), the area between the Mw ~ 9 large sliparea and M ~ 7 earthquake continued to creep even after
an M ~ 9 earthquake. In addition, our study showed agradual narrowing of the locked area during Stage II(Fig. 8), whereas Barbot (2020) reported continual creep
as in Stage I. This discrepancy is likely due to differences
in model settings between this study and Barbot (2020).
Barbot (2020) modeled the lower arc between the forearc
basement and the mantle wedge as a permanent structural boundary with a velocity-strengthening condition.
The assumption of a velocity-strengthening boundary
allows M ~ 7 interplate earthquakes to frequently occur
and always creep during interseismic periods, even aftera giant rupture. In addition, the locked region within theM ~ 9 source area did not change during an interseismic
period of an M ~ 9 earthquake cycle. That is, the physicalprocesses that shorten the time interval of MYG earthquakes differed between our results showing three stages
and those of Barbot (2020) without Stage II. If the rate at
which the locked area shrunk is faster than the simulation, a large afterslip would have occurred at the updip
area following the earthquake at the end of Stage I. The
locked area would then largely disappear. In that case,Stage I would be followed by Stage III. Seismic and geodetic observations on the updip side are important toknow if Stage II is present.
4.3 Multiple patch model with other M ~ 7 earthquakes

In the previous section, for simplicity, a circular unstable patch was located only in the Miyagi-ken-Oki segment. In reality, however, M ~ 7 earthquakes have been
observed on the updip side and, in some cases, occurred
after the Miyagi-ken-Oki earthquake (e.g., the earthquakes in 2005 and 2011, 1978 and 1981, and the 1930s
and 1939). We also conducted numerical simulations
using other models, which included seismic patches
that assumed the same radius at the same location as
SHL1 and SHL2 in the middle segment used in Nakata
et al. (2016). One of these is shown in Additional file 2:
Fig. S2.6 and Table S2. This model is the same as that inFig. 2 and Table 1, except for the SHL1 and SHL2 patches.
We performed numerical simulations for a periodof ~ 2700 years. The results showed a similar trend tothat in the model with only the MYG patch concerning
the spatiotemporal distribution of the M ~ 9 and MYG
area was also obtained in the model with two additional
patches. In the shallow part, the time development forthe slip deficit showed four M ~ 9 earthquakes (M = 8.97,

Page 13 of 16

8.95, 8.91, and 8.97) with 638, 553, and 643 years between
them, respectively (Additional file 2: Fig. S2.7). Earthquakes of M = 7.26–7.50 frequently occurred in the MYG
patch with time intervals of 50–136 years. As reported in
Sect. 3, even with this model, the time evolution of theslip deficit in the center of the MYG patch (Additional
file 2: Fig. S2.7) was not constant throughout an M9earthquake cycle. The slip deficit rate was large after theM ~ 9 earthquake but became low and constant as thenext M ~ 9 earthquake approached. The intervals between
each great earthquake and the subsequent M > 7 earthquake in the MYG patch were 28–56 years; these were
the shortest intervals during each M9 earthquake cycle.
Furthermore, the next MYG earthquake had the longest
interval from the previous one (> 100 years). At the latest
stage of the M9 earthquake cycle, the average recurrence
intervals of the four M > 7 earthquakes in the MYG patch
were 71–74 years. At the end of the M9 earthquake cycle,
the interval between the MYG and great earthquakes was
0.7–5.5 years in this model. The time interval between the
M = 7.43 earthquake in the MYG patch and the second
M ~ 9 earthquake (T = 1295 years, Cycle 2) was 5.5 years,
consistent with the observations (5.5 years between theM7.1 earthquake in 2005 and M9.0 mainshock).
At the SHL1 patch, two to three earthquakes of
M = 6.81–7.43 occurred during an M9 earthquake
cycle. At the SHL2 patch, zero to two M < 7 earthquakes
occurred during an M9 earthquake cycle. The time interval between the M = 7.43 earthquake in the SHL1 patch
and the first M ~ 9 earthquake (T = 657 years, Cycle 1)
was 55 h, consistent with observations from March 2011(51 h between the M7.3 foreshock and M9.0 mainshock).
The time interval and magnitude of the earthquakes that
occurred in the MYG and SHL1 patches in the final stageof the M9 earthquake cycle differed between cycles.
The time interval of each event just prior to an M ~ 9earthquake did not affect the behavior after the M ~ 9earthquake.
Characteristics for the slip deficit rate (Additional
file 2: Fig. S2.7), coseismic/afterslip distributions of the
M > 7 MYG earthquake (Additional file 2: Fig. S2.8), andcumulative slip distributions (Additional file 2: Fig. S2.9)
were similar to those in the one-patch model. A notable difference was that the aseismic slip (afterslip of an
MYG earthquake) that propagated to the east sometimes triggered another M > 7 and M ~ 7 earthquake at
the SHL1 and/or SHL2 patches within the M ~ 9 source,
respectively (Additional file 2: Fig. S2.8). In the present
simulations, there was also a case of an SHL1 earthquake
triggering an MYG earthquake. Thus, it is reasonable todiscuss characteristics common to all cycles concerning
MYG earthquakes based on the model without SHL1 andSHL2 patches.

Nakata et al. Progress in Earth and Planetary Science(2023) 10:34

In addition, we showed other models, which included
seismic patches that assumed the same radius at the same
location as FKS and IBK in the southern segment usedin the study by Nakata et al. (2016). Frictional parameters and temporal distribution of slip deficits are shown
in Additional file 2: Figs. S2.10, S2.11, and Table S3. Thismodel is the same as that in Fig. 2 and Table 1, except for
the SHL1, SHL2, FKS, and IBK patches. Spatial and temporal development of the slip deficit at the MYG patch
showed similar characteristics as the simple model.
4.4 March and May 2021 earthquakes

Next, we examined whether the 2021 M < 7 earthquakesthat occurred in the earliest stage of the great earthquake
cycle caused the transition from Stage I to II, or whether
Stage I continued. No M > 7 earthquakes have occurredin the region since 2011. As the source area was approximated using a single and homogeneous circle, we did notmodel M < 7 earthquakes similar to those in 2021 in the
MYG area. The 2021 earthquakes can simply be regarded
as small earthquakes that cannot be reproduced during
Stage I by this model.
The epicenters of the 2021 earthquakes were locatedat almost the exact site of the epicenter of the 1937
earthquake (Umino et al. 2007). In the later stage of theM ~ 9 earthquake cycle ending with the 2011 TohokuOki earthquake, the source area migrated from updip todowndip in the order of 1933 (Mj7.1), 1936 (Mj7.4), and1937 (Mj7.1) (Umino et al. 2007). In contrast, after the2021 earthquakes, seismic activity between the Marchand May source areas suggested that the afterslip of the2011 great earthquake continued in this area (Yoshida
et al. 2022). Thus, the western part of the 1978 source
does not appear to have fully locked after the 2021 earthquakes. Otherwise, the eastern part of the 1978 source
possibility remained locked. These situations differedfrom the beginning of Stage II in our simulation in thatthe whole MYG patch is locked. In addition, the interplate earthquakes indicated that the aseismic slip at the
downdip side of the MYG patch was still larger than
before 2011 (Fig. 7a; Yoshida et al. 2022).
Therefore, the current situation is still considered to bein the early-to-mid phase of Stage I. This suggests thatthe remaining areas of the 1978 source that had not ruptured in 2021 might rupture at the end of the remaining
short Stage I. Observing the spatial and temporal changes
in aseismic slips after the 2021 earthquakes in the eastern region are necessary to investigate the possibilityof future earthquakes in the eastern region of the 1978
earthquake rupture area (also those in the 1930s). It isalso important to determine whether the Tohoku-Oki
afterslip on the downdip side has continued. For more
realistic simulation results, we will adopt a multiscale

Page 14 of 16

circular patch model (Ide and Aochi 2013) to reproduce
a partial rupture in the MYG region, such as in the 1930s,
2005, and 2021 earthquakes.
4.5 Possibility of an M > 7 earthquake off Miyagi
within the next 30 years

The slips and average recurrence intervals of the simulated MYG earthquakes were larger and longer thanthose observed. As discussed in Sect. 4.2, we concluded
that the spatial relationship of slip and/or locked area
between the M9 earthquake and the MYG earthquakes
is significant in constraining the length of the recurrence
intervals. We next sought to examine the possibility of an
M > 7 earthquake occurring off Miyagi within the next
30 years based on our simulations. We would not showprobability (%) as HERP has presented because we do
not consider providing a percentage is the only direct
way to evaluate the long-term activity. Similar to Nakata
et al. (2016), the present results showed that “the timeinterval between the M ~ 9 earthquake and subsequent
M > 7 MYG earthquake was shorter than the average
recurrence time interval during the later stage of the M9earthquake cycle.” This might have been caused by stressaccumulation between the locked area at the updip sideand the large afterslip at the downdip side within andaround the MYG patch.
If a large amount of afterslip following an M ~ 9 earthquake occurred surrounding the MYG patch, it could
create a similar situation to that in Stage III. In that case,the time interval of the M > 7 MYG earthquake would be
relatively constant throughout the M ~ 9 cycle. However,this hypothesis was rejected given that the aseismic slip
observed after the 2011 Tohoku-Oki earthquake was sufficiently smaller on the updip side of the MYG patch thanthat on the downdip side (Fig. 7a).
Nakata et al. (2016) also presented several scenarios
in which the interval between an M ~ 9 earthquake anda subsequent M > 7 MYG earthquake is longer than the
average interval of several MYG earthquakes during thelater stage of the M9 earthquake cycle. A large slow slip
occurred in those scenarios instead of an M > 7 earthquake. The 2021 earthquakes differed from a simulated
slow slip in terms of slip velocity. Thus, scenarios with alonger time interval than the average recurrence interval after an M ~ 9 earthquake would be rejected. Hence,scenarios with a shorter time interval than the average
recurrence interval are still possible. Given that the current situation is considered to be in Stage I, an earthquake larger than those in 2021 will likely occur shortly.
These simulations showed that the variation in recurrence intervals within Stage III is small compared to thechanges from Stage I to Stage II and Stage II to Stage
III. Therefore, it can be concluded that the recurrenceNakata et al. Progress in Earth and Planetary Science(2023) 10:34

interval of MYG earthquakes in the later stage of an M9earthquake cycle is stable. Therefore, in the later stage,
it is reasonable to rely on the history of the several pastearthquakes regarding long-term evaluation, as HERP
had done prior to 2011. In contrast, in the earliest stageof the M9 earthquake cycle, long-term evaluation basedon the recent history is unsuitable because of exceptionalsituations based on the occurrence of the great earthquake. In this case, the long-term (30 years) evaluation
should refer to the spatiotemporal behavior of aseismic
slips in and around the MYG source in addition to thepast recurrence intervals.

5 Conclusions
Using numerical simulations, we investigated the spatiotemporal distribution of interplate slip velocity for the
M > 7 Miyagi-ken-Oki earthquake during an M ~ 9 earthquake cycle along the Japan Trench. Despite approximating the M > 7 source area with a single circle, thesimulation results were validated in terms of observed
characteristics of relatively constant time intervals for
the M > 7 Miyagi-ken-Oki earthquakes before the 2011
Tohoku-Oki earthquake, afterslip distributions of the
M > 7 Miyagi-ken-Oki earthquakes in 1978 and 2005,
coseismic slip distributions for the 2021 earthquakes
located in the western region of the primary rupture areasof the earthquakes in 1978, and repeating-earthquakederived interplate slip observations before and after the2011 Tohoku-Oki earthquake. These observations werequalitatively consistent with the simulations before andafter M ~ 9 earthquakes. Moreover, we reported that therecurrence intervals of the Miyagi-ken-Oki earthquakes
varied in the early and later stages of the M ~ 9 earthquake cycle. These differences might be related to theextent of the locked and afterslip areas in, and around,the Miyagi-ken-Oki source, particularly on the updip
side. Although the simulation results were not exactly thesame as the observations with respect to the time intervalof earthquakes and aseismic slip amounts, based on our
simulations, the recurrence interval of the Miyagi-kenOki earthquake was unstable during the early stage after
the M ~ 9 earthquake. For this period, observations ofaseismic slip and seismicity in and around the expectedsources are important to assess the timing of expected
Miyagi-ken-Oki earthquakes and long-term probability
evaluations based on past recurrence intervals.
Abbreviations
M	Magnitude
Mj	Magnitude as determined by the Japan Meteorological Agency
JMA	Japan Meteorological Agency
MYG	Miyagi-ken-Oki
HERP	Headquarters for Earthquake Research Promotion
BPT	Brownian Passage Time

Page 15 of 16

Supplementary Information
The online version contains supplementary material available at https://​doi.​
org/​10.​1186/​s40645-​023-​00566-y.
Additional file 1. 1.1-1.3. Detailed discription of Methods. Table S1.
Parameters used in this study. Figure S1. Schematic illustration of thedefinition of subfaults.
Additional file 2. Additional Figures and Tables. Figure S2.1. Positionof the squares used in Figure 7. Figures S2.2-2.5. Time developmentof simulated slip velocity during Cycle 3. Figure S2.6. Distribution offrictional parameters for three patches. Figure S2.7. Temporal distributionof slip deficits obtained using the friction model shown in Figure S2.6. Fig‑
ure S2.8. Time development of simulated slip velocity obtained using thefriction model shown in Figure S2.6. Figure S2.9. Cumulative slip in the
middle segment before and after an M ~ 9 earthquake obtained using thefriction model shown in Figure S2.6. Figure S2.10. Distribution of frictionalparameters for five patches. Figure S2.11. Temporal distribution of slip
deficits obtained using the friction model shown in Figure S2.10. Table S2.
Parameters for the northern, middle, southern, and deep segments andthe MYG, SHL1, and SHL2 patches shown in Figure S2.6. Table S3. Parameters for the FKS and IBK patches shown in Figure S2.10.
Acknowledgements
The computational resources of the Earth Simulator provided by JAMSTEC
and the Cyberscience Center, Tohoku University, were used for all numericalsimulations. Figures were drawn using Generic Mapping Tools (Wessel andSmith 1998). We would like to thank T. Matsuzawa, an associated editor, andthree reviewers for their thoughtful suggestions and insights. We would also
like to thank Editage (www.​edita​ge.​com) for English language editing.
Author contributions
RN performed the numerical simulations, proposed the topic, and conceived
and designed the study. NU helped in the preparation and interpretation of
the observation data. TH assisted with the preparation and interpretation of
the simulation results and collaborated with the corresponding author in theconstruction of the manuscript. RH helped validate the simulation results
based on the observational data and collaborated with the corresponding
author in the construction of the manuscript. All authors read and approved
the final manuscript.
Funding
This work was supported by JSPS KAKENHI (Grant Numbers JP21K04604,
JP19H05596, JP21K03713, and JP21H05206).
Availability of data and material
The datasets generated and/or analyzed during the current study are notpublicly available because of the large size of the datasets; however, the data
are available from the corresponding author on reasonable request.

Declarations
Competing interests
The authors declare that they have no competing interests.
Author details
1
Graduate School of Science, The University of Tokyo, 7‑3‑1 Hongo, Bunkyo‑Ku,
Tokyo 113‑0033, Japan. 2 Graduate School of Science, Tohoku University, 6‑6,
Aramaki‑Aza‑Aoba, Aoba‑Ku, Sendai 980‑8578, Japan. 3 Earthquake ResearchInstitute, The University of Tokyo, 1‑1‑1 Yayoi, Bunkyo‑Ku, Tokyo 113‑0032,
Japan. 4 Research and Development Center for Earthquake and Tsunami
Forecasting (FEAT), Research Institute for Marine Geodynamics (IMG), Japan
Agency for Marine-Earth Science and Technology (JAMSTEC), 3173‑25
Showa‑Machi, Kanazawa‑Ku, Yokohama 236‑0001, Japan.
Received: 29 September 2022 Accepted: 29 June 2023

Nakata et al. Progress in Earth and Planetary Science(2023) 10:34

References
Baba T, Ito A, Kaneda Y, Hayakawa T, Furumura T (2006) 3-D seismic wave velocity structures in the Nankai and Japan Trench subduction zones derived
from marine seismic surveys. Abstract S111–006 presented at the JpGU
Meeting 2006, Makuhari, Japan, 14–18 May 2006. https://​www2.​jpgu.​org/​
meeti​ng/​2006/​pdf/​S111/​S111-​006_e.​pdf. Accessed 1 July 2022.
Barbot S (2020) Frictional and structural controls of seismic super-cycles at
the Japan trench. Earth Planets Space 72:63. https://​doi.​org/​10.​1186/​
s40623-​020-​01185-3
Bird P (2003) An updated digital model of plate boundaries. Geochem Geophys Geosyst 4:1027. https://​doi.​org/​10.​1029/​2001G​C0002​52
Comninou M, Dundurs J (1975) The angular dislocation in a half space. J Elast
5:203–216. https://​doi.​org/​10.​1007/​BF001​26985
Dieterich JH (1979) Modeling of rock friction: 1. Experimental results andconstitutive equations. J Geophys Res 84:2161–2168. https://​doi.​org/​10.​1029/​JB084​iB05p​02161
Headquarters for Earthquake Research Promotion (2011) Evaluations of occurrence potentials of subduction-zone earthquakes. http://​www.​jishin.​go.​jp/​evalu​ation/​long_​term_​evalu​ation/​subdu​ction_​fault/. (Japanese).
Accessed 1 July 2022.
Headquarters for Earthquake Research Promotion (2021a) Evaluation of earthquake offshore Miyagi Prefecture on March 20, 2021a. https://​www.​jishin.​go.​jp/​main/​chousa/​21apr_​miyagi/​index-e.​htm. Accessed 1 July 2022.
Headquarters for Earthquake Research Promotion (2021b) Evaluation ofseismic activities for May 2021b. https://​www.​jishin.​go.​jp/​main/​chousa/​
21jun/​index-e.​htm. Accessed 1 July 2022.
Headquarters for Earthquake Research Promotion (2022) Evaluation of occurrence potentials of subduction-zone earthquakes. https://​www.​jishin.​go.​jp/​main/​chouk​ihyoka/​ichir​an.​pdf. Accessed 1 July 2022.
Hino R, Suzuki K, Yamamoto Y, Nishino M, Kanazawa T, Yamada T, Nakahigashi
K, Mochizuki K, Shinohara M, Kuwano A, Aoki G, Tanaka M, Araki E, KodairaS, Fujie G, Kaneda Y (2007) Aftershock distribution of the 2005 Miyagi-oki
earthquake (M 7.2) revealed by ocean bottom seismological observation:
preliminary results. Zisin 2(59):297–308. https://​doi.​org/​10.​4294/​zisin.​59.​
297
Hyodo M, Hori T, Kaneda Y (2016) A possible scenario for earlier occurrenceof the next Nankai earthquake due to triggering by an earthquake at
Hyuga-nada, off southwest Japan. Earth Planets Space 68(6):1–17. https://​doi.​org/​10.​1186/​s40623-​016-​0384-6
Ide S, Aochi H (2013) Historical seismicity and dynamic rupture process of the2011 Tohoku-Oki earthquake. Tectonophysics 600:1–13. https://​doi.​org/​
10.​1016/j.​tecto.​2012.​10.​018
Igarashi T, Matsuzawa T, Umino N, Hasegawa A (2001) Spatial distribution of
focal mechanisms for interplate and intraplate earthquakes associatedwith the subducting Pacific plate beneath the northeastern Japan arc: a
triple-planed deep seismic zone. J Geophys Res 106:2177–2191. https://​doi.​org/​10.​1029/​2000J​B9003​86
Iinuma T, Hino R, Kido M, Inazu D, Osada Y, Ito Y, Ohzono M, Tsushima H, Suzuki
S, Fujimoto H, Miura S (2012) Coseismic slip distribution of the 2011
off the Pacific Coast of Tohoku Earthquake (M9.0) refined by means ofseafloor geodetic data. J Geophys Res. https://​doi.​org/​10.​1029/​2012J​
B0091​86
Iinuma T, Hino R, Uchida N, Nakamura W, Kido M, Osada Y, Miura S (2016)
Seafloor observations indicate spatial separation of coseismic andpostseismic slips in the 2011 Tohoku earthquake. Nat Commun 7:13506.
https://​doi.​org/​10.​1038/​ncomm​s13506
Miura S, Iinuma T, Yui S, Sato T, Tachibana K, Hasegawa A (2007) Interplate slip
distribution inferred from co- and post-seismic deformation associatedwith the 2005 Miyagi-oki earthquake (N7.2). Zisin 2(59):371–379. https://​doi.​org/​10.​4294/​zisin.​59.​371
Matthews MV, Ellsworth WL, Reasenberg PA (2002) A Brownian model forrecurrent earthquakes. Bull Seismol Soc Am 92:2233–2250. https://​doi.​
org/​10.​1785/​01200​10267
Nakata R, Hori T, Hyodo M, Ariyoshi K (2016) Possible scenarios for occurrenceof M ~ 7 interplate earthquakes prior to and following the 2011 TohokuOki earthquake based on numerical simulation. Sci Rep 6:25704. https://​doi.​org/​10.​1038/​srep2​5704
Nakata R, Hori T, Miura S, Hino R (2021) Presence of interplate channel layer
controls of slip during and after the 2011 Tohoku-Oki earthquake
through the frictional characteristics. Sci Rep 11:6480. https://​doi.​org/​10.​1038/​s41598-​021-​86020-9

Page 16 of 16

Nakatani M (2001) Conceptual and physical clarification of rate and state friction: Frictional sliding as a thermally activated rheology. J Geophys Res
106:13347–13380. https://​doi.​org/​10.​1029/​2000J​B9004​53
Ozawa S, Nishimura T, Munekane H, Suito H, Kobayashi T, Tobita M, Imakiire
T (2012) Preceding, coseismic, and postseismic slips of the 2011 Tohoku
earthquake Japan. J Geophys Ress. https://​doi.​org/​10.​1029/​2011J​B0091​
20
Rice JR (1993) Spatio-temporal complexity of slip on a fault. J Geophys Res
98:9885–9907. https://​doi.​org/​10.​1029/​93JB0​0191
Ruina A (1983) Slip instability and state variable friction laws. J Geophys Res
88:10359–10370. https://​doi.​org/​10.​1029/​JB088​iB12p​10359
Suito H, Nishimura T, Tobita M, Imakiire T, Ozawa S (2011) Interplate fault slip
along the Japan Trench before the occurrence of the 2011 off the Pacificcoast of Tohoku Earthquake as inferred from GPS data. Earth Planets
Space 63:615–619
Uchida N, Bürgmann R (2021) A decade of lessons learned from the 2011
Tohoku-Oki earthquake. Rev Geophys. https://​doi.​org/​10.​1029/​2020R​
G0007​13
Uchida N, Matsuzawa T (2013) Pre- and postseismic slow slip surrounding the2011 Tohoku-oki earthquake rupture. Earth Planet Sci Lett 374:81–91.
https://​doi.​org/​10.​1016/j.​epsl.​2013.​05.​021
Uchida N, Matsuzawa T, Hirahara S, Hasegawa A (2006) Small repeating earthquakes and interplate creep around the 2005 Miyagi-oki earthquake (M= 7.2). Earth Planets Space 58:1577–1580
Uchida N, Matsuzawa T, Miura S, Hirahara S, Hasegawa A (2007) Interplate
quasi-static slip off Miyagi and Fukushima Prefectures estimated fromsmall repeating earthquake data. Zisin 2(59):287–295
Umino N, Kono T, Okada T, Nakajima J, Matsuzawa T, Uchida N, Hasegawa A,
Tamura Y, Aoki G (2007) Relocation of the M ~ 7 Miyagi-oki earthquakes
in the 1930s: Seismic slips of asperities that were re-ruptured during the1978 M7.4 Miyagi-oki earthquake? Zisin 2(59):325–337. https://​doi.​org/​10.​4294/​zisin.​59.​325
Utsu T (1979) Seismicity of Japan from 185 through 1925-a New Catalog
of Earthquakes of M > 6 Felt in Japan and smaller earthquakes whichcaused damage in Japan-. Bull Earthq Res Inst 54:253–308
Wessel P, Smith WHF (1998) New, improved version of generic mapping tools
released. Eos Trans Am Geophys Union 79:579. https://​doi.​org/​10.​1029/​
98EO0​0426
Yaginuma T, Okada T, Hasegawa A, Kato K, Takemura M, Yagi Y (2007) Coseismic slip distribution of the 2005 Miyagi-Oki (M 7.2) earthquake estimatedby inversion of strong-motion and teleseismic waveforms—Its relation
with the 1978 Miyagi-Oki earthquake (M 7.4). Zisin 2(60):43–53. https://​doi.​org/​10.​4294/​zisin.​60.​43
Yamanaka Y (2005) EIC Seismological_Note. No.168. https://​wwweic.​eri.u-​
tokyo.​ac.​jp/​sanchu/​Seismo_​Note/​2005/​EIC168.​html. Accessed 1 July
2022.
Yamanaka Y, Kikuchi M (2004) Asperity map along the subduction zone innortheastern Japan inferred from regional seismic data. J Geophys Res.
https://​doi.​org/​10.​1029/​2003J​B0026​83
Yoshida K, Matsuzawa T, Uchida N (2022) The 2021 Mw7.0 and Mw6.7 MiyagiOki earthquakes nucleated in a deep seismic/aseismic transition zone:
possible effects of transient instability due to the 2011 Tohoku earthquake. J Geophys Res 127:202. https://​doi.​org/​10.​1029/​2022J​B0248​87

Publisher’s NoteSpringer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.