Simulation of great earthquakes along the Nankai Trough: reproduction of event history, slip areas of the Showa Tonankai and Nankai earthquakes, heterogeneous slip-deficit rates, and long-term slow slip events – Earth, Planets and Space

In this section, we present and discuss the simulation results for the Nankai Trough plate boundary based on the best parameter distributions (Fig. 4; Table 2). To avoid the effect of initial values, we show only events after the first 600 years of the 5000-year simulation. We add the prefix “S” to elapsed time from the start of the simulation to avoid confusion with calendar dates. During S600–S5000, our simulation produced 106 seismic events with M_w between 6 and ~ 8.5 and repeating LSSEs (Additional file 1: Figs. S2 and S3). Among seismic events, 66 were great earthquakes or mega-earthquakes (M_w ~ 8 or higher), 31 were earthquakes of M_w ~ 7.5 in northern Hyuganada, and 9 were earthquakes of M_w 6-class in LSSE areas.

Earthquakes

Representative slip distributions

We used rupture patterns to classify our simulated 66 earthquakes of M_w ~ 8 or higher into the following six types (Fig. 5). Nine of them events could not be classified as one of these:

1. 1.

   The “Showa type” is a pair of earthquakes with epicenters off Cape Shiono in which the deeper parts of segments A and B rupture in an M_w 8.0–8.1 event within 3 years after segments C and D rupture in an M_w 7.9–8.1 event.

2. 2.

   The “Ansei type” is a pair of earthquakes in which segments A and B rupture in an M_w 8.3–8.4 event within 0.1 year after segments C–F rupture in an M_w 8.3–8.4 event.

3. 3.

   The “Hoei type” is the simultaneous (within 0.03 years) rupture of segments A–E in an M_w 8.4–8.6 event and segment Z in an M_w ~ 7.5 event.

4. 4.

   The “quasi-Showa type” is the same as the Showa type but includes a rupture in the shallower part of segment A.

5. 5.

   The “quasi-Hoei type” is the same as the Hoei type but includes a rupture of segments Z–E within 0.03–1 years.

6. 6.

   The “isolated Tonankai type” is an M_w 8-class Tonankai event (segments C and D) without an M_w 8-class Nankai event (segments A and B) within several tens of years.

Examples of slip distribution in six types of earthquake events in our 4400-year simulation. Stars indicate rupture initiation points. In each plot, the event number is at upper left (see Additional file 1: Fig. S2 for all 106 numbered events), and the numbers at lower right are the elapsed years in the simulation (S600–S5000), the moment magnitude, and the maximum slip displacement. Colors of the borders correspond to the event type

Following are some considerations that affected our definitions. Although the historic Ansei Tokai and Nankai earthquakes were about 32 h apart, in the definition of Ansei-type earthquakes we allowed the time lag to be as long as 0.1 year. With respect to Showa-type events, Baba and Cummins (2005) used tsunami data to estimate that the main rupture area of the 1946 Showa Nankai earthquakes was near the coast, and we targeted that area (target phenomenon 1e); however, the joint inversion analysis by Murotani et al. (2015) using seismic waveforms and geodetic data showed that an area of large slip on the Nankai segment amounting to M_w 8.4 was farther offshore. Accordingly, we defined an event with rupture area including offshore with magnitude less than M_w 8.4 as a quasi-Showa-type earthquake. With respect to Hoei-type events, the 1707 Hoei earthquake rupture may not have extended eastward to segment F (Matsu’ura 2012); thus, it may be reasonable to add the possibility that segments E and F remain unruptured to the definition of a Hoei-type earthquake given the low resolution of rupture areas inferred from historical documents. We therefore defined a rupture of segments Z–D or Z–E of magnitude M_w ≥ 8.5 within 1 year as a quasi-Hoei-type earthquake, although the recorded time lag between the Hoei Tokai/Tonankai and Nankai earthquakes is extremely short (Imai et al. 2011).

Figure 6 shows the simulated earthquakes during S2000–S3400, and Additional file 1: Figs S2 and S3 show events for the whole simulation period (S600–S5000). Our model showed that events east and west of Cape Shiono tend to occur in pairs, sometimes simultaneously but usually with time lags between them. The individual rupture areas of these paired events generally showed considerable variation (target phenomenon 1a), but even when two event pairs had identical rupture areas, they differed in their maximum displacements and moment magnitudes (compare events 23–24 and 32–33 in Additional file 1: Fig. S2).

Occurrence times and source regions of great Nankai Trough earthquakes during years S2000–S3400 of our simulation. Stars indicate rupture initiation points. Blue lines indicate rupture segments with large slips (> 5 m) near the trough axis. Red horizontal lines in area Z indicate Hyuganada earthquakes with M_w ~ 7.5; green lines in areas E and Z indicate M_w 6-class earthquakes. Numbers in parentheses indicate years between successive earthquakes. Annotation on the right side lists the magnitudes of Nankai and Tokai/Tonankai earthquakes; the superscripts indicate the event number, corresponding to the plots in Figs. 5 and Additional file 1: S2. Mega-earthquake magnitudes are shown in white on red background. The numbers in brackets indicate the time lag in years between the Tokai/Tonankai and Nankai earthquakes of event pairs; event pairs with superscripts b, c, and d mean time lags of 2 days, 1 day, and 2 weeks, respectively. The arrows indicate the occurrence order of great earthquakes, and the open arrows indicate the cases where the Nankai earthquake preceded. The earthquake types, shown on the right, are described in the text. Their colors correspond to Fig. 5. Question marks indicate rupture patterns that cannot be classified into any of six representative earthquake types shown in Fig. 5

Great earthquakes

Our model simulated six Showa-type earthquakes and reproduced the epicenters and rupture areas of the 1944 and 1946 Showa earthquakes as estimated from tsunami data (Baba and Cummins 2005) (compare the Showa type in Figs. 5 and 2a). There were eight quasi-Showa-type earthquakes. There were six Ansei-type earthquakes, and our model roughly reproduced the coseismic slip area of the 1854 Ansei earthquake pair (Aida 1981a, b) (compare Figs. 5, 2b). There were two Hoei-type earthquakes, and our model roughly reproduced the coseismic slip area of the 1707 Hoei earthquake (Furumura et al. 2011; Matsu’ura 2012; Kobayashi et al. 2018) (compare Figs. 5 and 2c). There were three quasi-Hoei-type earthquakes.

The sequence of events 95–103 in S4572–S4799 (Additional file 1: Figs. S2 and S3) was generally consistent with observations in terms of the event order and slip areas of the historical (Hoei, Ansei, and Showa) earthquake sequences, if the quasi-Hoei type is included. In this sequence, great earthquakes had occurrence intervals of 137 and 86 years, similar to the historical intervals of ~ 150 and 90 years (target phenomenon 1b). The interval between events of the Showa-type pairs (events 102 and 103) in the simulation was 2.6 years, which is consistent with the historical interval of ~ 2 years (target phenomenon 1c). Although the interval between events in the simulated Ansei-type pairs (events 99 and 100) differed from the historical record by an order of magnitude (0.1 y in the simulation vs. 32 h), the simulation is qualitatively consistent with the historical record in that the difference was shorter for Ansei-type than for Showa-type event pairs (target phenomenon 1c). However, the interval between the quasi-Hoei-type event pair (events 95 and 96) was 0.7 years, which is longer than the 0.1-year interval of the Ansei-type pairs.

The 66 great earthquakes included 31 pairs in which events east and west of Cape Shiono were separated by 4.5 years or less. Among these pairs, the eastern event occurred first in 28 pairs and the western event occurred first in 3 pairs (events 12 and 13, events 40 and 41, and events 44 and 45 in Additional file 1: Fig. A). In addition, three events had simultaneous ruptures of both the east and west sides of Cape Shiono (events 10, 15, and 48 in Additional file 1: Fig. S2). Each of these three events initiated west of Cape Shiono, and their magnitudes were M_w 8.5–8.6. One event, an M_w 8.1 isolated Tonankai earthquake (event 39 in Additional file 1: Fig. S2) may be regarded as similar to the 1498 Meio Tokai earthquake, because the existence of its western counterpart is doubtful (e.g., Ishibashi 1998).

Earthquake Research Committee (2013) evaluated the occurrence potentials of Nankai Trough earthquakes with the use of a time-predictable recurrence model (Shimazaki and Nakata 1980). Hashimoto (2022) reexamined the long-term evaluation of the occurrence of great earthquakes along the Nankai Trough issued by Earthquake Research Committee (2013) and related materials, and found several problems in the adoption of the time-predictable model (e.g., neglect of measurement errors of uplift during historical earthquakes; inconsistency between assumed uplift and geodetic or geomorphological data). Thus he recommended to adopt a simple average of recurrence intervals instead of the time-predictable model to calculate the probability of the next earthquake. Hori et al. (2012) showed a time-predictable behavior in a simple hierarchical asperity model consisting of two small patches, but also pointed out that the behavior did not necessarily appear in general. Figure 7 shows the cumulative slip displacements at points 1–5 in Fig. 4f. The Tonankai asperity (area including point 3), where ruptures tended to initiate first, fit neither the time-predictable nor slip-predictable models. Therefore, it appears that the recurrence interval changed in response to local differences in stress conditions in surrounding asperities. The slip-predictable model was a better fit for the other four asperities (although the point 1 asperity was rather periodic), which tended to slip in association with ruptures of adjacent asperities. From these results, when discussing the probability of occurrence of the next great earthquake, it is important to evaluate without persisting in the time-predictable model.

Time evolution of cumulative slip displacement (thick line) at locations marked by points 1–5 in Fig. 4f. Numbers in parentheses are event numbers, as shown in Figs. 6 and Additional file 1: S2. Upper and lower thin lines indicate slip-predictable and time-predictable models, respectively

Mega-earthquakes

The 1707 Hoei earthquake, at M_u 8.6, is generally considered to be the largest known earthquake in the Nankai Trough (Table 1). Our simulation produced eight mega-earthquakes (events 2, 10, 15, 26, 36, 48, 88, and 96 in Additional file 1: Fig. S2), with occurrence intervals of 287–1654 years, and the average interval of 560 years is similar to the recurrence intervals of 400–600 years (target phenomenon 1d) estimated from field investigations (Maemoku 1988; Shishikura et al. 2008; Okamura and Matsuoka 2012). The slip was large in the area of the subducting seamount southeast of Cape Muroto with the low slip-deficit rate, even when we assigned the area a large (L) value to produce stable sliding.

We also examined events with large slips (> 5 m) near the trough axis (blue lines in Figs. 6 and Additional file 1: S3). About 40% (= 4/11) of the Tokai/Tonankai events that ruptured segment F reached the trough axis. Slips were large at the trough axis when the area with large ({sigma }^{mathrm{eff}}) (point 5 in Fig. 4e) and the area with the subducting seamount (Tosabae) both ruptured (e.g., event 88 in Additional file 1: Fig. S2; compare events 85, 88, 90; event 85 did not rupture both areas and event 90 ruptured point 5 but not the Tosabae seamount).

Our findings also offer insight into the role of deep locked zones. Okamura and Shishikura (2020) proposed, based on geological data on crustal movements of the forearc wedge, that the 1854 Ansei Tokai earthquake was ruptured with the deep locked zone at ~ 20 km depth beneath the Akaishi Mountains in the Tokai region, whereas the 1944 Showa Tonankai and 1946 Showa Nankai earthquakes did not rupture the deep locked zone. With respect to this proposal, we were able to reproduce the variations in the eastern end of the ruptures of historic (1707, 1854, and 1944) earthquakes by setting barrier-like asperities with larger ({sigma }^{mathrm{eff}}) and ({h}^{*}) at 5–20 km depth in the area including points 1 and 2 in the Tokai district (Fig. 4e, f). They also proposed that mega-earthquakes with a recurrence interval of 400–600 years, such as the 1707 Hoei earthquake, are controlled by the locked zone at 20–30 km depth beneath the Kii Mountains in Kii Peninsula. However, our simulation reproduced mega-earthquakes by setting barrier-like asperities in the area including point 5 off Shikoku. The placement of a strong locked zone beneath the Kii Mountains is an open question.

Time lag between Tokai/Tonankai and Nankai earthquakes

The Working Group to Investigate Disaster Prevention Based on Monitoring and Assessment of Seismicity along the Nankai Trough (2017) considered the possibility of pairs of great earthquakes occurring in close succession on opposite sides of Cape Shiono. The working group found that among the 96 events of M_w 8.0 or more that occurred worldwide after 1900, 10 were followed within 3 days by an earthquake of the same magnitude in an adjacent region. They showed that this fraction of event pairs could be approximated by the Omori–Utsu equation (Utsu 1961).

Our simulation produced 34 pairs of Tokai/Tonankai and Nankai great earthquakes on adjacent fault segments, of which 3 were simultaneous and 30 involved a time lag of less than 3 years (1 pair, events 30 and 31, had a time lag > 3 year). The time lags tended to cluster around peaks at several days and 0.1, 0.4, 0.7, and 2.6 years (Additional file 1: Fig. S3), a result that is inconsistent with the Omori–Utsu equation (Utsu 1961). Our model may simply differ from that equation in this respect, or there may be time lags that are inherent in the Nankai Trough region. Although the sample size is very small, the four historical earthquake pairs on the Nankai Trough had lags of either a few days (1–3 August 1361 for the Shohei par and 23–24 December 1854 for the Ansei pair) or a few years (17 December 1096 to 22 February 1099 for the Kowa/Eicho pair and 7 December 1944 to 21 December 1946 for the Showa pair).

Hyuganada earthquakes

Our model simulated M_w 7.4–7.6 events in the northern part of the Hyuganada segment with recurrence intervals of 81–292 years (Figs. 6; Additional file 1: S3). The time interval between the two Hyuganada events in the historical record (Fig. 1), ~ 260 years, lies within this range (target phenomenon 2). Aseismic slips, which were also simulated in the asperity area in northern Hyuganada (point 6 in Additional file 1: Fig. S1), would release stress in that area.

Hyodo et al. (2016) showed that whether a great earthquake was triggered by a Hyuganada earthquake depended on the timing of the Hyuganada earthquake. However, because they set a shorter distance between the asperities of Nankai and Hyuganada earthquakes than we did in our model (their Nankai earthquake asperity reached Cape Ashizuri, unlike ours in Fig. 4f), the two asperities in their model could easily influence each other. Although not showing in a figure, when we ran a simulation with a shorter distance between our two asperities (areas including points 4 and 6 in Fig. 4f), it produced frequently Nankai earthquakes that also ruptured the Hyuganada asperity, a result inconsistent with the historical record. Under the set of spatial parameters necessary to reproduce the historical record, shown in Fig. 4f, no great earthquake initiated at the asperity in northern Hyuganada. Although Hyodo et al. (2016) also suggested that the 1498 Meio Nankai earthquake might have been triggered by a Meio Hyuganada earthquake, both of these earthquakes are doubtful (Ishibashi 1998; Harada et al. 2017). Our results indicate that it is difficult for Hyuganada earthquakes to trigger the Nankai earthquakes. Note that the seismic gap with a recurrence interval of ~ 260 years may be an apparent one due to the paucity of historical records (see “Slip-deficit rate distribution“).

M
_w 6-class earthquakes in LSSE areas

In our model, nine M_w 6-class earthquakes ruptured the deeper part of the seismogenic zone in segments E and Z (events 17, 21, 28, 42, 43, 46, 50, 63, and 98 in Additional file 1: Fig. S2; green lines in Figs. 6 and Additional file 1: S3). According to the seismic catalog, which includes estimated magnitudes for pre-instrumental events (Usami (2003) for events before 1885, Utsu (1982) for 1885–1922 events, and Japan Meteorological Agency (2017) for 1923–present events), 7 earthquakes of M 6-class occurred in the Tokai district and the Bungo Channel during 1600–2015 (Fig. 8). Although the hypocenter locations of events before 1884 have limited precision, these events might have ruptured the deeper part of the seismogenic zone, similar to the events simulated by our model.

Epicentral distribution and magnitude diagram of M 6-class earthquakes in the Tokai district (solid stars) and the Bungo Channel (open stars) during 1600–2015. D Indicates the depth of an earthquake

In our simulation, M_w 6-class earthquakes or acceleration of slip in the Bungo Channel LSSE area sometimes induced M_w ~ 7.5 Hyuganada earthquakes. For example (see Additional file 1: Fig. S2), M_w 6-class earthquakes in the Bungo Channel LSSE area appeared only seven times in 4400 years (events 17, 28, 43, 46, 50, 63, and 98), and for four pairs (events 17–18, 28–29, 46–47, and 50–51) of them, the time differences between the M_w 6-class earthquakes and the M_w ~ 7.5 Hyuganada earthquakes were within 0.04 years. In addition, simulated M_w ~ 7.5 Hyuganada earthquakes occurred 31 times in 4400 years, and five (events 22, 54, 66, 77, and 94) of them were earthquakes whose rupture starting point (shown by a star in Additional file 1: Fig. S2) was located in the Bungo Channel LSSE area. In the Tokai LSSE area, in contrast, the slip here did not propagate to the adjacent area and M_w 6-class earthquakes (events 21 and 42 in Additional file 1: Fig. S2) in the Tokai LSSE area did not induce any larger event. This is because we modeled subducting ridges by assigning a large ({sigma }^{mathrm{eff}}) (Sect. “Tuning of effective normal stress ({sigma }^{mathrm{eff}})“) and large (L) (Sect. “Tuning of characteristic displacement L“) to the seismogenic zone (areas including points 1 and 2 in Fig. 4) to suppress the rupture propagation during the 1944 Tonankai earthquake. This modeling also behaved as a barrier against slip in the Tokai LSSE area to propagate into the subducting ridge area, resulting in M_w 6-class earthquakes (events 21 and 42) staying only in the Tokai LSSE area.

In our model, the ruptures of great earthquakes tended to initiate off Cape Shiono (Additional file 1: Figs. S2 and S3). This tendency was an effect not only of the plate configuration (Takayama et al. 2008), but also of the setting of a large (L) off Cape Shiono in order to reproduce the rupture starting points of the 1944/1946 Showa earthquakes (Additional file 1: Section S2.1) and low slip-deficit rates (Additional file 1: Section S2.4). Because the Tokai and Bungo Channel LSSE areas are far from the initiation points of great earthquakes, there may be no cases of M_w 6-class earthquakes in either of these LSSE areas triggering great earthquakes.

Slip-deficit rate distribution

Offshore observations have recently revealed that slip-deficit rates are heterogeneous on the Nankai Trough plate boundary (Yokota et al. 2016; Nishimura et al. 2018). In this section, we compare slip-deficit rate distributions between our model and observations (Nishimura et al. 2018) (Fig. 2d) by using the sequence of events during S4338–S4491, which is similar to the historical sequence of the 1854 Ansei, 1944/1946 Showa, and 1968 Hyuganada earthquakes. Nishimura et al. (2018) used inland GNSS data from 2005 to 2009 and seafloor GNSS/A data from 2004 to 2016 in a block modeling approach to estimate slip-deficit rates on the plate boundary by superimposing rigid motions between blocks and elastic deformations due to slip defects between blocks. Because more than 70 years have passed since the 1946 Showa Nankai earthquake, we examined the spatial distributions of the slip-deficit rate on the plate interface at ~ 70, 100, and 130 years after the Nankai earthquake (event 93 at S4428) in our simulation (Fig. 9). The slip-deficit rate is the difference between the plate convergence rate (Fig. 4c) and the slip rate in each cell of the model. In the northern Hyuganada segment, an M_w 7.5 event occurred in S4491, 63 years after event 93, indicating that afterslip around the Hyuganada focal region had temporarily decreased the slip-deficit rate to nearly zero off Cape Ashizuri (Fig. 9a). Then the slip-deficit rate continued to decrease in the Hyuganada focal region while it recovered slightly in the surrounding area (Fig. 9b, c). Note that the time difference between the Nankai and Hyuganada events in our model (63 years) was longer than that in the historical record (22 years). Although the slip-deficit rate distributions in this example from our model were generally larger than those observed and larger slip-deficit rates extended further up-dip (< 10 km depth) than the observation, their spatial patterns were generally similar at all three time points and relatively consistent to the observation.

Slip-deficit rates on the plate interface at a ~ 70, b 100, and c 130 years after the Nankai earthquake in S4428 (event 93) in our simulation. Solid and broken black lines indicate the 3.0 and 2.4 cm/y contours (standard error < 2.0 cm/y), respectively, of currently observed slip-deficit rates (Nishimura et al. 2018). Numbers at lower right indicate the elapsed time from the start of the simulation, the elapsed time from the previous earthquake, and the leading time to the next earthquake

The 4400-year simulation produced several examples of Showa-type event (red outlines in Additional file 1: Fig. S2) that were followed by a great earthquake pair, four that ruptured segments A–D (quasi-Showa type) and one that ruptured segments Z–D (quasi-Hoei type). In these cases, an eastern event of M_w 8.0 and a western event of M_w 8.3 or more occurred with a time lag of about half a year. Figure 10 shows the slip-deficit rate distributions and slip velocity distributions before and after the quasi-Hoei-type earthquake event (events 95–97) that followed a Showa-type event (events 92 and 93). The slip-deficit rates before event 95 in S4572, which ruptured the Tonankai segment, were decreasing slightly, then changed significantly near the starting point, from 3 cm/y to nearly zero off Cape Shiono between S4553 (Fig. 10b) and S4571 (Fig. 10d). After event 95, the interplate coupling in the Tonankai area and the Kii Channel LSSE area recovered, while slip near the starting point of event 96 in S4573, which ruptured the Nankai segment, progressed up-dip and nucleated (see the velocity distribution). The evolution of interplate coupling before and after the quasi-Showa-type events that followed a Showa-type event was similar (for example, events 75–76 and 78–79). The current slip-deficit rate off Cape Shiono is estimated to be almost zero (Nishimura et al. 2018) (Fig. 2d). The high slip-deficit rate region in the Tonankai segment is currently smaller (Fig. 2d) than that in the simulation result for S4572 at the time of event 95 (Fig. 10f), and the next great event has not yet happened. This result suggests that the actual nucleation area may be larger than the area used in our model, with a critical nucleation size of about 27 km (point 3 in Table 2).

Distributions of (left) slip velocity and (right) slip-deficit rate on the plate interface southeast of Cape Shiono at six points in time before and during the quasi-Hoei-type earthquake (events 95–97) which followed the Showa-type earthquake of events 92 and 93. Slip velocity is normalized by the plate convergence rate. The quasi-Hoei-type event (see Fig. 5) consisted of a Tonankai event (95), Nankai event (96) and Hyuganada event (97). Plots depict a, b ~ 20 years before and c, d ~ 1 year before event 95; e, f the start of event 95 (epicenter shown by star) g, h immediately and i, j ~ 270 days after event 95, and k, l the start of event 96 (epicenter shown by star)

Although Nishimura et al. (2018) reported an area of high slip-deficit rate beneath the Bungo Channel, that rate might decrease in the long term because no M_w 7-class LSSE was included in the observation data. When we assigned a large ({sigma }^{mathrm{eff}}) to northern Hyuganada (point 6 in Fig. 4e) to reproduce the seismic gap of ~ 260 years, the simulated slip-deficit rate was higher than the observations. Because this difference in slip-deficit rates was large, the apparent absence of earthquakes in the observations may reflect the paucity of the historical record.

Other studies besides Nishimura et al. (2018) have estimated slip-deficit rates on the Nankai Trough plate interface (Noda et al. 2018; Kimura et al. 2019). They reported slip-deficit rate distributions along the strike direction similar to that of Nishimura et al. (2018). Kimura et al. (2019) also estimated the slip-deficit rate from land-based GNSS and seafloor GNSS/A data using a block motion model similar to that of Nishimura et al. (2018), adding observation data from five seafloor GNSS/A stations near Suruga Bay, Kumano Basin, and the trough axis managed by Nagoya University. To make the spatial uncertainty of their estimated values uniform, they included subfaults with low spatial resolution beneath the ocean, reflecting the low spatial density of the seafloor GNSS/A network. This made it possible to estimate slip-deficits near the trough axis, albeit at a lower spatial resolution. High slip-deficit rates were estimated to the south of Cape Ashizuri, the Kii Channel, and the Shima Peninsula. However, these regions have been the site of shallow LSSEs (Yokota and Ishikawa 2020) and VLFEs (Takemura et al. 2019), indicating low interplate coupling. In an inversion analysis incorporating viscoelastic response to estimate slip-deficit rates from land-based GNSS and seafloor GNSS/A data, Noda et al. (2018) reported a peak in the slip-deficit rate that was closer to the trough axis than that of Nishimura et al. (2018). Whereas Nishimura et al. (2018) and Kimura et al. (2019) used the same plate configuration model as our study, Noda et al. (2018) used a CAMP model (Hashimoto et al. 2004) that was based mainly on intraplate earthquakes in the Philippine Sea slab. Watanabe et al. (2018) noted that differences in the plate configuration can affect estimates of slip-deficit rates. In sum, we consider the results of Nishimura et al. (2018) to be suitable for comparison with our simulation results.

We confirmed that large (L) has a certain effect in reproducing low slip-deficit rate distributions, but actual observation results show a much smaller slip-deficit rate southeast of Cape Muroto. To better reproduce the observations, we could set (left(a-bright){sigma }^{mathrm{eff}}) close to zero by defining local areas with extremely small ({sigma }^{mathrm{eff}}) or by making (left(a-bright)) positive. Shallow VLFEs (Takemura et al. 2019) (Fig. 2d), which indicate the presence of fluid, occur where seamounts and ridges have subducted, indicating locally small ({sigma }^{mathrm{eff}}) (Fig. 2d). Hori (2006) and Kodaira et al. (2006) assigned (a-b>0) to a wedge-shaped area off Cape Shiono to represent the unlocked plate interface based on the existence of the highly fractured oceanic crust caused by strike–slip faults indicating stress release through gradual deformation there. It might be appropriate to do so in other low slip-deficit rate regions if strike–slip faults suggesting the unlocked plate interface are present.

LSSEs

Figure 11 shows four time series of simulated slip velocity in LSSE areas, in which intermediate (unsaturated) peaks in slip velocity correspond to LSSEs. Beneath the Tokai district (point 8 in Fig. 4e), simulated LSSEs occurred in the latter half of interseismic periods. Beneath the Kii Channel (point 9), the simulated LSSE amplitude was very small, being large only after the unpaired Tonankai earthquake in S2169 (event 39). Beneath western Shikoku (point 10), LSSEs tended to be excited by M_w ~ 7.5 earthquakes in northern Hyuganada or M_w ~ 6 earthquakes in the Bungo Channel LSSE area. Observations based on GNSS data (Kobayashi 2017; Takagi et al. 2019) indicate the occurrence of western Shikoku LSSEs within one year after Bungo Channel LSSEs. Our simulation appears to reproduce the response of western Shikoku LSSEs to surrounding stress disturbances (although the cause of disturbance is different such as LSSEs in observations and earthquakes in our model). In the Bungo Channel district (point 11), simulated LSSEs were nearly constant throughout the interseismic periods.

Time evolution of slip velocity in the a Tokai, b Kii Channel, c western Shikoku, and d Bungo Channel segments (points 8–11 in Fig. 4e, f). Unsaturated low-velocity fluctuations represent LSSEs. Saturated high velocities (> 25 cm/y) represent seismic slips induced by M_w 8-class great earthquakes and M_w ~ 7.5 Hyuganada earthquakes in the surrounding area, and there are only nine M_w 6-class earthquakes with spontaneous seismic slips in the Tokai and Bungo LSSE areas (Additional file 1: Fig. S2). Dashed lines indicate the plate convergence rate at each point. Labels R (for rupture) in a and d indicate great earthquakes in the seismogenic zone up-dip from the LSSE areas. Blue and orange dashed lines between c and d indicate M_w ~ 7.5 earthquakes in area Z and M_w 6-class earthquakes in the Bungo LSSE area, respectively. Detailed data for the shaded time periods are presented in Fig. 12

Beneath the Tokai and Bungo Channel districts, simulated LSSEs tended not to occur after ruptures up-dip in the seismogenic zone (R in Fig. 11). The stress heterogeneities associated with a non-ruptured seismogenic zone cause LSSEs to be generated in the deeper zone, as pointed out by Hirose and Maeda (2013). In fact, in the Tokai region (segments E and F), the area up-dip of the Tokai LSSE area did not rupture in the 1944 Tonankai earthquake. Likewise, in the Bungo Channel region (segment Z), because the 1946 Nankai earthquake rupture did not extend past Cape Ashizuri, the area up-dip of the Bungo Channel LSSE area did not rupture.

Figure 12 shows detailed time series of simulated LSSEs in the 100-year periods shown shaded in Fig. 11. The recurrence interval of the simulated Tokai LSSEs was 11.3–11.8 years, and their amplitudes increased with time until the M_w 6-class earthquake (event 21) in S1484. The recurrence interval of the simulated Kii Channel LSSEs was 8.8–14.4 years, and their amplitudes varied. Prior to event 40, the slip accelerated at the lower end of the Nankai asperity (blue line in Fig. 12b), the pre-slip progressed to the shallow plate interface, and the rupture started near the trough axis (Additional file 1: Fig. S2). The increase in LSSE amplitude several years before event 40 was due to the pre-slip at the lower end of the asperity, and it is unlikely that it directly triggered event 40. The recurrence interval of the simulated western Shikoku LSSEs was 2.9–7.7 years, and their intervals lengthened and their amplitudes decreased with time. The simulated Bungo Channel LSSEs differed before and after the S4711 Hyuganada earthquake (event 101). Before the earthquake, the variations of the recurrence interval (6.0–6.1 years) and the amplitude were small, and afterward the amplitude first became large and then decreased with time, whereas the occurrence interval was initially short and then increased with time.

Details of Fig. 11 showing the time evolution of slip velocity during 100-year periods of the simulation in the a Tokai, b Kii Channel, c western Shikoku, and d Bungo Channel segments (shaded periods in Fig. 11). Earthquakes (EQ) are labeled with red lines and event numbers are in parentheses. Circles denote events illustrated in Fig. 13. Blue line in b indicates the time evolution of slip velocity at point X in Fig. 13b

Because Bungo Channel LSSEs before 1980 could not be estimated, even with pre-GNSS geodetic data (Kobayashi and Yamamoto 2011), we could not confirm velocity changes of the western Shikoku and Bungo Channel LSSEs before and after the 1968 Hyuganada earthquake. In general, the occurrence intervals of the simulated LSSEs were consistent with the observations (Ozawa 2017). Note that the simulated LSSEs did not directly trigger great Nankai earthquakes, because they were too distant from the earthquake epicenters.

The critical nucleation size ({h}^{*}) in the Tokai LSSE area (point 8) is slightly smaller than those in the other LSSE areas (Table 2). Since the ratio of the asperity size to the critical nucleation size is large, unstable sliding is likely to occur (Kato 2003). As a result, Tokai LSSEs may have tended to increase gradually in amplitude, and in some cases seismic slip may have occurred (Fig. 12a). The recurrence interval of earthquakes is roughly regulated by (left(b-aright){sigma }^{mathrm{eff}}) of the seismogenic zone (e.g., Stuart 1988) and the plate convergence rate (e.g., Hirose and Maeda 2013). Because the plate convergence rate is slower in the east than in the west (Nishimura et al. 2018) (Fig. 4c), recurrence intervals should be longer in the east if the friction parameters are the same. Comparing the three segments other than Tokai, which have the same friction parameters (Fig. 12), it can be confirmed that the recurrence interval is slightly longer in the Kii Channel, in the east. The nucleation size also depends on the aspect ratio of a rectangular asperity, and when the ratio (alpha =w/l), where (w) and (l) indicate the fault width and length, is large, the nucleation size becomes small (Kato 2003). Although it is a qualitative consideration, because the shape of the LSSE area is not a simple rectangle, the Kii Channel LSSE area (point 9) should tend more toward unstable sliding than the Bungo Channel LSSE area (point 11), yet the seismic slip occurred in the Bungo Channel. Differences in the nucleation size of the LSSE area may not be the only defining characteristics of LSSEs, because all LSSEs tend to end with a transition to stable sliding at the plate convergence speed and because LSSEs are excited by ambient stress heterogeneity and disturbances. These differences also probably affect the LSSE velocity and recurrence interval in each area.

Figure 13 shows the spatial distributions of slips on the plate interface during the LSSEs marked with circles in Fig. 12. We estimated LSSE magnitudes within areas with small ({sigma }^{mathrm{eff}}) and (L) (Fig. 4e, f) by summing over the period when the slip rate exceeds the plate convergence rate at the reference points for those areas (points 8–11). The magnitudes of the simulated LSSEs in the Kii Channel and western Shikoku segments were consistent with the observations (see Fig. 6 in Ozawa 2017), but those in the Tokai and Bungo Channel segments were underestimated compared to the observations. Hirose and Maeda (2013) showed that large values of (left(b-aright){sigma }^{mathrm{eff}}) or small (L) values lead to LSSEs with large magnitudes. However, large values of (left(b-aright){sigma }^{mathrm{eff}}) are inappropriate because they also lengthen the recurrence interval of LSSEs and are inconsistent with the observation. In our model, although LSSE amplitude increased as (L) decreased, the calculations tended to become unstable as (L) approached the size of the computation cells (see “Computational instability due to cell size“). Therefore, (L) cannot be reduced any further given the resolution of our model. Because the calculation cost increases with the resolution of the model, reproducing the magnitudes of LSSEs in these areas must await future studies. Ochi and Kato (2013), who did not eliminate steady-state fluctuations associated with plate subduction, have suggested that the moment magnitude of the 2000–2005 Tokai LSSE was 6.6, not 7.2; in that case, the simulated Tokai LSSE would be consistent with observations. Because we estimated LSSE magnitudes simply from the amount of slip during the period, it is appropriate to compare our results with those of Ochi and Kato (2013). Note that the magnitude of Bungo Channel LSSEs is the same with or without steady-state fluctuations (Ochi 2015). The durations of simulated LSSEs were about 1–2 years in all areas (Fig. 13), shorter than the observed durations (Fig. 3) except in the Bungo Channel segment.

Distribution of slip on the plate interface during LSSEs. Black lines denote the edges of the source area. LSSEs in a–c and d correspond to LSSEs indicated by circles in a–c and d in Fig. 12, respectively. Other symbols are as shown in Fig. 2

Figure 14 plots the evolving relationship between the slip velocity and the shear stress during various events in the four LSSE areas. In the Tokai LSSE area (Fig. 14a), when the area up-dip of the Tokai LSSEs did not rupture during event 19 (note that the curve begins just before event 19), the shear stress was almost constant as the velocity dropped sharply and the plate interface stuck firmly once. After that, LSSEs occasionally released the accumulated stress, as shown by the tight loops in the curve. On the other hand, after event 13, which ruptured the area up-dip from the LSSE area, the shear stress increased as the velocity decreased and then stagnated near the plate convergence rate, and stress was released through stable sliding, such that LSSEs did not occur.

Trajectory diagrams showing the relation between slip velocity and shear stress in the a Tokai, b Kii Channel, c western Shikoku, and d Bungo Channel segments (points 8–11 in Fig. 4e, f). Each line corresponds to the simulation period shown at upper left and includes the number of the event during that period, shown in Additional file 1: Fig. S2

In the Kii Channel segment (Fig. 14b), the behavior after the isolated Tonankai earthquake (event 39) was similar to the case in which LSSEs occurred in the Tokai LSSE area. Under the influence of the adjacent Nankai asperity (point 4) up-dip from the Kii Channel LSSE area, which did not rupture, the Kii Channel LSSE area stuck firmly once, and then stress was released by occasional LSSEs. On the other hand, in a case of sequential rupture of the Tonankai and Nankai segments (events 52 and 53), after the Tonankai earthquake (event 52), the shear stress remained almost constant as the velocity dropped sharply and the plate interface stuck firmly once, as in the case of event 39. However, the Nankai earthquake (event 53) soon followed, and then the velocity stagnated near the plate convergence rate and stress was released by stable sliding, such that LSSEs did not occur.

In the western Shikoku segment (Fig. 14c), after the Nankai earthquake (event 53), the velocity stagnated near the plate convergence rate and entered a stable sliding cycle (see Figs. 11, 12). This is the reason why the Nankai asperity (point 4) up-dip from the LSSE area ruptured during every Nankai earthquake (the timing of the large slip velocity > 25 cm/y in Fig. 11c is the same as the rupture of the Nankai asperity). However, LSSEs occurred as a result of the stress disturbance caused by an M_w ~ 7.5 earthquake in northern Hyuganada (blue lines in Fig. 11c or event 54 in Fig. 12c) or an M_w ~ 6 earthquake in the Bungo Channel LSSE area (orange dashed lines in Fig. 11c). Their amplitude tended to decrease with time.

In the Bungo Channel segment (Fig. 14d), as in the Tokai LSSE area, LSSEs occasionally occurred when the area up-dip was unruptured (event 103) and did not occur when it was ruptured (events 96 and 97). After an M_w ~ 6 earthquake in the Bungo Channel LSSE area (event 98), the behavior was the same as when the up-dip area was unruptured. That is, the potential for LSSEs was increased by the combination of an unruptured up-dip area, as in the Tokai LSSE area (Fig. 14a), and the stress disturbance from middle-class earthquakes, as in the western Shikoku LSSE area (Fig. 14c). Observations (Fig. 3) show that LSSE activity is greater in the Bungo Channel than in other LSSE areas. From the above, our simulation shows that LSSEs are excited by heterogeneity and disturbances in ambient stress, and that it is important that the plate interface sticks firmly once under the conditionally unstable state (where (a-b<0) and the LSSE area size is comparable to ({h}^{*})).

Among simulated LSSEs beneath the Kii Channel (point 9), only the set of LSSEs after the S2169 unpaired Tonankai earthquake was large, which was caused by a heterogeneous stress distribution due to an unruptured up-dip area (adjacent Nankai asperity) like Tokai and Bungo Channel LSSEs. However, in fact the 1946 Nankai earthquake ruptured the area up-dip from the LSSE area (Baba and Cummins 2005; Murotani et al. 2015) (Fig. 2a) and LSSEs have occurred beneath the Kii Channel at least since 1996 (e.g., Kobayashi 2014, 2017; Ozawa 2017). Accordingly, our model did not simulate them. Although we cannot clarify the reason for this at present, one possible explanation is that an unknown heterogeneous stress distribution was generated beneath the Kii Channel at the 1946 Nankai earthquake.

The configuration of the Philippine Sea slab varies greatly. Mitsui and Hirahara (2006), modeling the effect of slab dip angle on LSSE occurrence, pointed out that lateral changes of slab dip are important factors controlling LSSE occurrence. After that, however, the increase in knowledge of plate configuration and LSSEs reveals that LSSEs also occur in areas where dip changes in the strike direction of the slab are small, such as western Shikoku and Bungo Channel (Fig. 4e, f). Consequently, lateral changes of slab dip are a factor affecting LSSE distributions, but not a necessary condition for LSSEs.

With regard to the relationship between the slab configuration and fluid distributions, Morishige and van Keken (2017) showed in a 3-D fluid model that fluids released from the slab are localized depending on the slab configuration, and they suggested that fluids beneath the Kii Channel tend to migrate away. Seismic tomography results (see cross section in Fig. 6 of Hirose et al. 2008) indicate that there is less fluid beneath the Kii Channel than beneath the Tokai district and the Bungo Channel. In this study, we assigned the same effective normal stress ({sigma }^{mathrm{eff}}) to the Kii Channel, western Shikoku, and Bungo Channel LSSE areas, but it may be that the Kii Channel should have a larger ({sigma }^{mathrm{eff}}) value. However, since ({sigma }^{mathrm{eff}}) defines the LSSE recurrence interval (Hirose and Maeda 2013), it is constrained by the observed values (Fig. 3) to a limited range of variation.

Although LSSE areas occur in a consistent depth range, (left(b-aright){sigma }^{mathrm{eff}}) must change greatly along the trough axis because LSSE areas are intermittent along the axis. Because frictional parameters (a) and (b) depend on depth (Blanpied et al. 1998), (left(b-aright)) should vary little along the trough axis; accordingly, it would be appropriate to allow ({sigma }^{mathrm{eff}}) to be heterogeneous, as in this study. However, it has been pointed out that (left(a-bright)) also depends on ({sigma }^{mathrm{eff}}) (Sawai et al. 2016) and sliding velocity (Boulton et al. 2018). Therefore, both ({sigma }^{mathrm{eff}}) and (left(b-aright)) would change along the trough axis with variations in fluid presence, even if the depth of the LSSE zone is constant. To produce more realistic models, parameters will need to be chosen with this in mind.

Problems to address in future studies

It is well known that a non-planar fault induces changes in normal stress that can significantly alter the fault strength under both static and dynamic loading (e.g., Dieterich and Smith 2009). In addition, Linker and Dieterich (1992) described RSF with a different state evolution law for normal stress changes. In this study, the triangular cells were arranged to fit the configuration of the plate boundary, but we did not consider the temporal variation of normal stress for the following reasons. Takayama et al. (2008) showed that, owing to shear slip, the value of the slip response function is 1–3 orders of magnitude smaller in the normal direction than the shear direction for the same 3-D plate configuration used in this study. That is, the temporal variation of normal stress due to shear slip is small compared to that of shear stress. Moreover, incorporating the temporal variation of normal stress doubles the calculation cost, which is not realistic for parameter tuning by trial and error, as we did in this study. Note that previous studies have not considered the temporal variation of the normal stress in simulations that introduced 3-D plate configurations in the Nankai Trough (Hirose and Maeda 2013; Hyodo and Hori 2013; Nakata et al. 2014; Hyodo et al. 2016). This remains an issue for future work.

The quasi-dynamic simulation used in this study approximates the stress change associated with the radiation of seismic waves by the effect of the radiation damping term (the second term on the right side of Eq. (1)). On the other hand, fully dynamic simulations fully incorporate stress changes associated with the propagation of seismic waves (e.g., Lapusta et al. 2000; Kaneko et al. 2010). Since fully dynamic simulation is computationally very expensive (Thomas et al. 2014), it is not realistic to apply it to a 3-D model of the entire Nankai Trough. Thomas et al. (2014) reported that fully and quasi-dynamic simulations generated similar slip patterns, although the slip velocity, the rupture propagation speed, and the degree of rupture stopping in the velocity-strengthening region during an earthquake were quantitatively different. Lapusta and Liu (2009) showed that a quasi-dynamic simulation set to the seismic radiation correction factor (eta) = 4 in Eq. (1) had high similarity to a fully dynamic simulation, but differed from it in detail with respect to the slip velocity, rupture propagation speed, and slip amount during an earthquake; they also found that the rupture pattern in the quasi-dynamic simulation appeared monotonous compared to that in the fully dynamic simulation. In this study, we set (eta) = 1 according to previous studies (e.g., Hori 2006; Hyodo and Hori 2013; Hirose and Maeda 2013; Nakata et al. 2014; Hyodo et al. 2016). When rescaled with (eta), the behaviors of (eta) = 1 and (eta) = 4 during an earthquake match well (Lapusta and Liu 2009). Therefore, even if (eta) = 4 is set in this study, the rupture pattern will not change significantly qualitatively. Note that a quasi-dynamic simulation with (eta) = 3 was performed off Tohoku by Nakata et al. (2016).

There is growing evidence that friction may be much lower at seismic slip velocities than RSF laws predict (e.g., Di Toro et al. 2004; Rice 2006). When this seismic weakening mechanism was introduced into the friction constitutive law, the behavior was significantly different between fully and quasi-dynamic simulations, and the rupture pattern became more complicated in the latter (Thomas et al. 2014). In the case of a simple two-dimensional model where two velocity-weakening regions were embedded in the velocity-strengthening region, earthquakes that ruptured both of them at once were repeated in a fully dynamic simulation, whereas small earthquakes that ruptured part of one velocity-weakening region were frequent in the quasi-dynamic simulation while large earthquakes that ruptured the entire area were rare. Accordingly, if we had incorporated a seismic weakening mechanism in our quasi-dynamic simulation, even if we had assigned the same parameters to the regions including points 1–3 and points 4–5 in Fig. 4f, we might have been able to reproduce the complex rupture patterns shown in Fig. 1. On the other hand, because a fully dynamic model with a seismic weakening mechanism produces a very simplified rupture pattern, parameter distributions must be more complicated than what we used to reproduce these rupture patterns. This is a future issue.

We targeted rupture areas (Fig. 1) on the basis of the Earthquake Research Committee (2013) compilation, but different interpretations are possible. Seno (2012) proposed a model in which Ansei- and Hoei-type earthquakes alternate, Ansei-type earthquakes rupturing segments A, B, and D–F (Fig. 1) at intervals of 356–412 years, and Hoei-type earthquakes rupturing segments A–D at intervals of 237–474 years. According to that model, the next great Nankai Trough earthquake would not occur for another 200 years. The Hoei-type earthquake of Seno (2012) corresponds to our Showa-type and quasi-Showa-type earthquakes, which were common in our simulation; however, our simulation did not produce any of Seno’s Ansei-type earthquakes because it was not one of our targeted phenomena. To do so, our model might have to allow effective normal stress to vary with time, as Seno (2012) suggested.

In our model, the eastern edges of Tokai/Tonankai earthquake ruptures did not vary much after S3000 (Additional file 1: Fig. S3), a result suggesting that this pattern may constitute a limit cycle. In addition, except for the Nankai earthquakes of events 88 and 96, the sequence of quasi-Showa-, Ansei-, and Showa-type events occurred repeatedly. This condition might be fundamental to actual great earthquakes along the Nankai Trough. By adding perturbations representing phenomena that we did not incorporate previously, our model might reproduce, more complicated and realistic phenomena. For instance, SSSEs of M_w 5–6 have also occurred at the edge of the seismogenic zone (Nishimura et al. 2013), and the interaction with great earthquakes is a concern. An M_w 5.8 earthquake on 1 April 2016 appeared to occur on the plate boundary off Cape Shiono. Because this location is near the starting point of the 1944 Showa Tonankai earthquake, an M_w 6-class asperity should be included in the future simulations with our model. In addition, Pollitz and Sacks (1995) estimated that the M_w 8.0 1891 Nobi earthquake in the continental plate beneath central Japan has delayed the anticipated Tokai earthquake by at least 20 years. Incorporating stress perturbations associated with surrounding earthquakes (e.g., Kuroki et al. 2004) would make our simulation model more realistic.

As discussed in “Plate configuration” section, we had to modify the near-surface domain of our model because the upper end of the plate boundary in the Nankai Trough axis is at a water depth of 4–5 km. In an alternative model incorporating the effect of surface topography, the large calculation errors in very shallow regions are partially obviated by assigning positive (a-b) values, which produce steady slip, to those regions (Ohtani and Hirahara 2015). However, analysis of a core sample collected off Cape Shiono has shown that coseismic slip propagated to the trough axis at least once (Sakaguchi et al. 2011). Therefore, it is necessary to assign negative (a-b) to the trough axis, as we did, and errors in shallow cells would therefore affect the slip behavior. Although setting a smaller cell size would improve the model accuracy, it would also increase computational costs. Consequently, our model did not incorporate surface topography.

It has been pointed out that the configuration of the eastern end of the Philippine Sea plate, east of 138°E and north of 35°N in the Tokai region, is shallower than that used in this study (Matsubara et al. 2021). Although differences in plate configuration are known to affect simulation results (e.g., Takayama et al. 2008), the relatively low plate convergence rate in this region (Fig. 4c) means the effect on our results should be small.

We modeled general aspects of historical earthquakes, including seismic and aseismic slips, to some extent by applying constant (a) and (a-b) and heterogeneous ({sigma }^{mathrm{eff}}) and (L) to the model space, and presented the best model for reproducing the historical events. Different results can be obtained by using other parameter distributions (see Additional file 1: Section S2). For example, when we changed (L) from the best model value of 0.22 m to 0.25 m (a change in critical nucleation size from 62.9 km to 71.4 km) at point 2 off the Tokai district, variation was suppressed at the eastern edge of simulated Tokai/Tonankai earthquake ruptures. As another example, even the small change of making the parameters of the LSSE areas beneath the Kii Channel and western Shikoku the same as those in the surrounding areas affected the simulated earthquake record. Thus, the precise simulation of recurring LSSEs beneath all of southwestern Japan, including Kyushu, is essential to improve models of great earthquakes. This in turn requires improvements in the monitoring technology for LSSEs.

In this study, we focused on reproducing various phenomena on the Nankai Trough plate boundary. However, most geodetic stations are on land, and it remains necessary to verify whether stations on land can capture changes of state at any location on the plate boundary. In addition, it is essential to construct a simulation model that can simulate various phenomena comprehensively, including crustal deformation, before it will be possible to attempt precise forecasts of the next great Nankai earthquake.

Although many offshore seismic stations have been emplaced in recent years, the spatial density of offshore stations is much lower than that of inland stations. A high-density network similar to DONET (Kaneda et al. 2015; Kawaguchi et al. 2015) on the seafloor of the Nankai Trough would be very useful for monitoring earthquakes, tsunamis, and crustal deformation in real time. Denser station coverage on the seafloor will improve the resolution of interplate coupling estimates, which will lead to more realistic models of the Nankai Trough plate boundary.