Stress change in southwest Japan due to the 1944–1946 Nankai megathrust rupture sequence based on a 3-D heterogeneous rheological model

Abstract

The Nankai trough has repeatedly experienced large megathrust earthquakes at intervals of 100–200 years. The inland region of southwest (SW) Japan has a seismically active period from 50 years before to 10 years after megathrust earthquakes. To assess the activities of inland earthquakes after megathrust earthquakes, we need to quantitatively evaluate the postseismic stress accumulation on the inland source faults considering plausible viscoelastic relaxation. Recent studies have shown the importance of low-viscosity layers along the lithosphere-asthenosphere boundary (LAB layer) in postseismic deformation. In the present study, we focus on the most recent ruptures, the 1944 M7.9 Tonankai and the 1946 M8.0 Nankai earthquakes, estimating the 4-year stress change on the source faults in SW Japan in a forward modeling approach. The 1944–1946 megathrust rupture sequence was followed by severe ~ M7 inland earthquakes, such as the 1945 M6.8 Mikawa and 1948 M7.1 Fukui earthquakes. For stress calculation, we used a highly detailed finite element model incorporating the actual topography and the plausible viscoelastic underground structure from past studies. The calculated inland stress field shows the dominance of the coseismic change during the 1944 and 1946 earthquakes and little contribution from viscoelastic relaxation. In contrast, viscoelastic relaxation has a significant effect on stress in the slab, indicating the importance of quantifying the viscosity of the LAB layer. Based on the calculated stress, we evaluated the change in the Coulomb failure stress (ΔCFS) on each source fault. The ΔCFS is generally positive on the strike-slip faults east of 135°E due to the 1944 rupture. In contrast, the ΔCFS on the faults west of 135°E, including the Median Tectonic Line segments, became positive due to the 1946 rupture. The occurrence of the damaging earthquakes in 1945 and 1948 can be explained by the calculated ΔCFS. The ΔCFS on the recent earthquake faults of the recent damaging earthquakes such as the 2016 M7.2 Kumamoto earthquake is generally negative, suggesting the delay in stress accumulation. The ΔCFS on the source faults of the intra-slab earthquakes differ significantly as large as tens of kilopascals depending on the viscosity of the LAB layer.

Graphical Abstract

Introduction

Under the southwest (SW) Japan arc, the Philippine Sea (PHS) plate has been subducting at the Nankai trough, causing repeating megathrust ruptures at intervals of 100–200 years (e.g., Ando 1975). According to historical documents, inland earthquakes in SW Japan (~ 200 km far from the trench) become active in the period from ~ 50 years before a megathrust event to ~ 10 years after an event (Utsu 1974; Shimazaki 1976; Seno 1979; Mogi 1981; Hori and Oike 1996). For quantitative assessment of such inland earthquake activity, we need to evaluate how stress accumulates in the inland source faults due to the cycle of megathrust ruptures. For the period before the coming megathrust rupture, we can calculate stress accumulation rates from a nearly steady megathrust locking pattern (slip deficit pattern) (e.g., Nishimura 2018; Saito et al. 2018). For the period after a megathrust rupture, however, the stress change strongly depends on the slip distribution and its viscoelastic effect. Thus, we need to evaluate several stress change patterns depending on possible rupture scenarios in advance.

The Nankai megathrust cycle has shown a variety of rupture patterns. The most recent is the 1944 and 1946 series of ruptures (the M7.9 Tonankai and M8.0 Nankai earthquakes), where the eastern and western halves separately ruptured at 2-year intervals. This rupture sequence was followed by damaging inland earthquakes of about M7, such as the 1945 M6.8 Mikawa and 1948 M7.1 Fukui earthquakes. In particular, the 1945 Mikawa earthquake occurred immediately after the first half rupture in 1944 and before the second rupture in 1946, indicating the importance of estimating the individual impacts of each rupture as well. Many studies have shown stress change due to the megathrust ruptures in the Nankai trough. Hori and Oike (1999) and Ogata (2004) showed the relation between the 1944–1946 rupture and the inland earthquake occurrence by stress calculation based on the elastic half-space model. Pollitz and Sacks (1997) computed the effect of the 1944–1946 rupture on the occurrence of the 1995 M7.3 Kobe earthquake by considering viscoelastic relaxation of the lower crust. Using a simplified elastic–viscoelastic layered half-space earth model, Shikakura et al. (2014) computed inland stress development in the central part of SW Japan considering more than 1000 years of megathrust rupture cycles along the Nankai trough. Their results showed the general occurrence of reverse-type earthquakes before the megathrust ruptures and strike-slip earthquakes after the megathrust ruptures. Mitogawa and Nishimura (2020) also computed the stress change on the source faults in SW Japan considering the periodic rupture of a large rectangular fault for the megathrust in an elastic–viscoelastic layered half-space.

Recently, an integrated onshore-offshore geodetic dataset from the 2011 M9 Tohoku earthquake, northeast Japan, revealed that postseismic displacement observed at the surface were strongly influenced by viscoelastic relaxation in the asthenosphere (Sun et al. 2014; Watanabe et al. 2014). Sun et al. (2014) and Freed et al. (2017) demonstrated the necessity of using a plausible viscoelastic structure in subduction zones, such as continental-oceanic viscosity contrasts and low-viscosity layers along the lithosphere-asthenosphere boundary (LAB layer) to reproduce the offshore postseismic movement. Using a two-dimensional model, Muto et al. (2013, 2016) showed the effect of the detailed heterogeneous viscosity structure of the arc on the postseismic displacements. Hashima et al. (2014) and Becker et al. (2018) demonstrated how postseismic stresses changed with viscoelastic relaxation due to the Tohoku earthquake. In South America, Li et al. (2018) showed the viscosity variation in the mantle under the continent affect the 8-year postseismic displacements due to the 2010 Maule earthquake. In the SW Japan arc, Suito (2017) suggested a viscosity structure similar to that of northeast Japan through his analysis of the 2004 M7.1 and M7.4 Kii Peninsula earthquakes. So far, no study has systematically shown possible postseismic loading effects on the intra-plate source faults due to contribution from the plausible viscosity structure stated above, particularly including the LAB layer.

Realistic stress calculations due to such a megathrust rupture require models with detailed surface topography and 3-D internal structure of large study area. However, such models have significant computational costs because they require very small elements to obtain reliable solutions. The development of analysis methods for crustal deformation, which effectively make use of recent high-performance computers, has enabled numerical computations of this kind (Fujita et al. 2017). Based on this method, Ichimura et al. (2016) and Yamaguchi et al. (2018) developed a high-performance finite element (FE) method for viscoelastic crustal deformation. Murakami et al. (2024) applied this high-performance model to calculate postseismic displacements due to a M9-class earthquake scenario under the Nankai trough and discussed the ability of the integrated onshore-offshore geodetic data to constrain viscosity structure, particularly the LAB layer.

In this study, we calculated the 4-year stress change in SW Japan due to the 1944–1946 rupture sequence in a forward modeling approach. We used the high-fidelity FE method by Ichimura et al. (2016) to evaluate the possible effect of viscoelastic relaxation in the viscoelastic regions under the Nankai trough, particularly the LAB layer, with high numerical accuracy. Then we computed the coseismic and postseismic change in the Coulomb failure stress (ΔCFS) on the 142 source faults in SW Japan. We directly compared the results of our stress predictions with the occurrence of postseismic damaging earthquakes, such as the 1945 Mikawa and 1948 Fukui earthquakes to validate our results. We also examined postseismic stress loading on recent earthquakes as the 1995 M7.3 Kobe and the 2016 M7.2 Kumamoto earthquakes to see the difference with postseismic damaging earthquakes.

Tectonic setting

Figure 1a shows the configuration of tectonic plates around SW Japan. Subduction of the PHS plate has formed the SW Japan and Ryukyu arcs. The SW Japan arc corresponds to the subduction of relatively younger part of the PHS plate (27–15 Ma, Okino et al. 1994). The direction of relative motion around SW Japan is oblique (N55°W), which is considered to drive the right-lateral strike slip of the Median Tectonic Line (MTL) (e.g., Sangawa 1977; Sato et al. 2015). The tectonic environment as described above continues to influence the fault development and present seismic activity in SW Japan.

Fig. 1

a Tectonic plates around the SW Japan arc. Solid lines with triangles are trenches. The blue line is the Median Tectonic Line (MTL). Moment tensors of Mw ≥ 6.5 and depth ≤ 30 km since 1976 are from the Global Centroid Moment Tensor Catalog (Dziewonski et al. 1981; Nettles et al. 2011; Ekström et al. 2012). Quadrangle shows the area of (b–d). b Rectangular fault model west of 138°E from HERP (2009b), F-net Catalog, Kakehi and Iwata (1992), Sagiya et al. (2002), and Suito (2017). Parameters of fault geometry and their correspondence with fault are shown in Additional file 1: Table S1 and Fig. S1. Faults shown by bold lines with year denote the source fault for the damaging earthquakes in Table 2, on which the 4-year stress change are shown in Figs. 9, 10, 11 and Additional file 1: Figs. S4–S6. c M ≥ 6.5 intra-plate earthquakes west of 138°E between 1854 and 1944 megathrust earthquakes. Color in symbols indicates elapsed year from 1854. d M ≥ 6.5 intra-plate earthquakes west of 138°E after 1944 megathrust earthquake. Color in symbols indicates elapsed year from 1944. Squares are megathrust earthquakes, circles are inland earthquakes, and triangles are intra-slab earthquakes. Earthquakes shown by thick lines are listed in Table 2. PAC, Pacific plate; EUR, Eurasia plate; PHS, Philippine Sea plate; SWJA, Southwest Japan arc; RA, Ryukyu arc; IBA, Izu-Bonin arc; NT, Nankai trough; RT, Ryukyu trench; KS, Kyushu; SB, Shikoku basin; OT, Okinawa trough; KP, Kii peninsula; TB, Tosa Bay

Figure 1b and Additional file 1: Table S1 show the distribution of the source faults in SW Japan, which have been active in the recent few hundred thousand years (HERP 2009b). Figure 1b clearly shows that source faults are mainly developed in the inner zone. Generally, the region east of ~ 135°E is east–west compressional due to the subduction of the Pacific plate and the collisional effect of the Izu-Bonin arc. Reverse and conjugate strike-slip faults corresponding to the stress field have developed. Between 132°E and 135°E, the MTL develops, as can be seen in the topography. Fewer source faults have been identified in this region than in the neighboring regions. West of 131°E is a transitional zone to the Ryukyu arc. Here, extensions of the volcanic front and the backarc basin of the Ryukyu arc reach Kyushu Island (Tada 1984, 1985), which causes increased volcanic activity. In particular, a prominent north–south extensional structure with numerous normal faults, called the Beppu-Shimabara Graben, develops in the middle part of Kyushu (Matsumoto 1979).

Under the Nankai trough, megathrust earthquakes occur at intervals of 100–200 years (e.g., Ando 1975). The rupture patterns of megathrusts are diverse; while the 1707 earthquake ruptured the whole locked region, the 1854 earthquake ruptured the eastern half, after which the western half ruptured 1 day later (Ando 1975; Ishibashi 2004), and the 1944–1946 rupture had a time lag of 2 years. The slip distribution of the 1944 Tonankai and 1946 Nankai earthquakes have been estimated by many studies with different methods and datasets, such as strong ground motion, tsunami, and surface deformation (e.g., Ando 1975, 1982; Yabuki and Matsu’ura 1992; Sagiya and Thatcher 1999; Tanioka and Satake 2001; Ichinose et al. 2003; Baba et al. 2002; Baba and Cummins 2005; Murotani et al. 2015; Sherrill and Johnson 2021). These results commonly show a slip distribution near the southeast offshore area of the Kii Peninsula for the 1944 rupture, and from Tosa Bay to the trench axis area on the south side of the Kii Peninsula in the 1946 rupture, though the details of these results differ depending on the characteristics of the analyzed data and the distribution of observation points.

As the 2011 M9 Tohoku earthquake clearly showed, megathrust earthquakes strongly affect the earthquake activity in neighboring regions by static and dynamic stress change (e.g., Asano et al. 2011; Hirose et al. 2011; Ishibe et al. 2011, 2015; Toda et al. 2011; Hasegawa et al. 2012; Fukushima et al. 2018; Zhang et al. 2023). In SW Japan, inland earthquakes are activated from ~ 50 years before to ~ 10 years after a megathrust earthquake (Utsu 1974; Shimazaki 1976; Seno 1979; Mogi 1981; Hori and Oike 1996). Figure 1c, d illustrate the occurrence of M ≥ 6.5 intra-plate earthquakes west of 138°E both in the overriding and subducting plates for the period between the 1854 and 1944 megathrust events and the period after the 1944 event, respectively, based on the catalog of Utsu (1990, 2002, 2004, The later updates are added by the International Institute of Seismology and Earthquake Engineering), F-net Moment Tensor Catalog to complement the earthquakes after 2013, and Headquarters for Earthquake Research Promotion (HERP 2009a). Hereafter, we refer to intra-plate earthquakes in the overriding and subducting plates as inland earthquakes and intra-slab earthquakes, respectively. In Fig. 1c, d, the color of the plotted earthquakes indicates the elapsed time from the last megathrust earthquakes in 1854 and 1944, respectively. Note that we assume the same recurrence time for the occurrence of the next megathrust rupture following Hori and Oike (1996). Thus, we use the same color scales for the elapsed time in both periods (Fig. 1c, d). We can notice the relative scarcity of the inland earthquakes in the period of 10–40 years both after the 1854 and 1944 megathrust earthquakes. We may call the earthquakes before and after the quiet period the postseismic and interseismic earthquakes, respectively. Figure 1c, d also show that the region of the postseismic earthquake occurrence are relatively limited around 136–137°E and 132°E and the interseismic earthquakes occurred more widely west of 136°E, suggesting a difference in postseismic and interseismic stress accumulation patterns. As for the intra-slab earthquakes, we cannot find a clear tendency between postseismic and interseismic patterns due to their scarcity.

Methods

Finite element model

We used the viscoelastic analysis technique using FE method developed by Ichimura et al. (2016), and the FE method improved by Hori et al. (2021) and Murakami et al. (2024). This study’s model and the analysis method are the same as those used by Murakami et al. (2024).

We assumed the Maxwellian viscoelasticity for our problem. The governing equations are written as follows:

$${\sigma }_{ij,j}+{f}_{i}=0,$$

(1)

$${\dot{\sigma }}_{ij}=\lambda {\dot{\varepsilon }}_{kk}{\delta }_{ij}+2G{\dot{\varepsilon }}_{ij}-\frac{G}{\eta }\left({\sigma }_{ij}-\frac{1}{3}{\sigma }_{kk}{\delta }_{ij}\right),$$

(2)

$${\varepsilon }_{ij}=\frac{1}{2}\left({u}_{i,j}+{u}_{j,i}\right).$$

(3)

Here, σ_ij, ε_ij, u_i, and f_i represent the stress, strain, displacement, and body force, respectively, with dot denoting the derivative in time. λ, G, and η are the first Lamé parameter, shear modulus, and viscosity, respectively. δ_ij is the Kronecker delta. The Einstein summation convention was used.

We take the origin of our model at the 0-m height point at (133°E, 33.5°N) on the ellipsoidal earth model GRS80 and the Cartesian x-, y-, and z-axes in the eastward, northward, and upward directions, respectively. Then, we take the model region in − 1248 km ≤ x ≤ 1248 km, − 1248 km ≤ y ≤ 1248 km, and − 1100 km ≤ z ≤ 0 km. The topography, internal elastic structure, and geometry of the plate boundary are based on the Japan Integrated Velocity Structure Model (JIVSM) (Koketsu et al. 2009, 2012).

Our model considers the viscosity of three regions: the mantle wedge, the LAB layer, and the oceanic mantle under the LAB layer (Fig. 2c). Note that the boundary between the mantle wedge and the oceanic mantle is taken along a horizontal plane extending from the tip of the slab at a depth of 75 km for convenience of model construction, which differs from the models in Sun et al. (2014) and Suito (2017). The overriding plate and subducting slab are assumed to be elastic. The elastic thickness is set as 45 km based on the maximum depth of the cold wedge in Sherrill and Johnson (2021). Also, the thickness of the slab is 30 km based on Sherrill and Johnson (2021). We assume that the LAB layer is 10 km thick. Following Sun et al. (2014) and Suito (2017), we assume five cases for the viscosity structure (Table 1): uniform viscosity (1 × 10¹⁹ Pa s, 2 × 10¹⁸ Pa s) (Cases 1 and 5), continental-oceanic viscosity contrast (2 × 10¹⁸ Pa s, 1 × 10¹⁹ Pa s) without the LAB layer (Case 2), continental-oceanic contrast with a moderately low LAB layer (2.5 × 10¹⁸ Pa s) (Case 3), and continental-oceanic contrast with an extremely low LAB layer (2.5 × 10¹⁷ Pa s) (Case 4).

Fig. 2

a Finite element model used in this study. b Close up view of the meshes around plate boundary. c Viscosity structure assumed in this study. The viscosity in each region is shown in Table 1

Table 1 Viscosity of the mantle wedge, low viscosity layer along the lithosphere-asthenosphere boundary (LAB layer), and oceanic mantle under the LAB layer used in the calculation

The analysis is performed using the FE method based on a tetrahedral quadratic element mesh to accurately consider the complex geometry and heterogeneous material properties of crustal structures as well as stress-free boundary conditions on the ground surface. The model region is automatically divided into unstructured meshes of quadratic tetrahedral elements (Ichimura et al. 2016). The minimum size of these elements is 2000 m except that the elements around the slip region are set as 500 m for numerical convergence (Fig. 2b). Then, our model consists of 1.0 × 10⁹ elements with 4.2 × 10⁹ degrees of freedom (three times the number of nodes). The boundary condition on the ground surface is free except including gravitational effects due to the vertical displacements (Barbot and Fialko 2010). The sides and bottom of the model are fixed to 0 under the Dirichlet boundary condition.

The coseismic slip of the 1944 and 1946 earthquakes are given with the slip node technique (Melosh and Raefsky 1981). To minimize the computation cases but also consider the variety of slip distributions seen in past studies, this study expressed the slip distribution U(x,t) by the superposition of five basic slip areas u_i (x) (i = 1, 2, 3, 4, 5) (Fig. 3) as

$$U\left( {{\bf{x}},t} \right) \, = \, \left[ {a_1 u_1 \left( {\bf{x}} \right) \, + a_2 u_2 \left( {\bf{x}} \right)} \right]H\left( {t-T_{1944} } \right) \, + \, \left[ {a_3 u_3 \left( {\bf{x}} \right) \, + a_4 u_4 \left( {\bf{x}} \right) \, + a_5 u_5 \left( {\bf{x}} \right)} \right]H\left( {t-T_{1946} } \right),$$

(4)

where u₁ (x) and u₂ (x) denote the 1944 slip and u₃ (x), u₄ (x), and u₅ (x) denote the 1946 slip with weighting factors a =  (a₁, a₂, a₃, a₄, a₅). H(t) is the Heaviside function, which takes a value of 0 for t < 0 and 1 for t > 0. T₁₉₄₄ = 1944.94 years and T₁₉₄₆ = 1946.97 years. Each u_i (x) (i = 1, 2, 3, 4, 5) takes a value of 1 m in the central part and decays sinusoidally outward (Fig. 3). The slip direction is N55°W following Miyazaki and Heki (2001) for all u_i (x). We computed the 4-year response to each basic slip u_i (x) with constant time steps of dt = 20 days to include the 1945 Mikawa and 1948 Fukui earthquakes in the period of stress computation. Following the latest study on the 1944–1946 slip distribution from Sherrill and Johnson (2021), we take a = (a₁, a₂, a₃, a₄, a₅) = (2, 5, 8, 0, 8).

Fig. 3

a Five basic slip regions. Sections AA′ and BB′ are sections under which stress fields are shown in Figs. 6, 7, Additional file 1: Figs. S2 and S3. b Schematic illustration of the slip distribution u_i(x) in each basic slip region

Under the above-stated setting, we solved the viscoelastic equation system (1–3) by the algorithm of Ichimura et al. (2016). Because our FE model has a large number of degrees of freedom, it is difficult to calculate the long-term viscoelastic response in a feasible amount of time using conventional computation methods. In this study, we performed the above-stated crustal deformation analysis in a feasible computation time by using a computation method that can efficiently utilize GPU computers, which have become widely used in recent years. Here, we use a highly efficient parallel iterative solver based on the conjugate gradient method with preconditioning based on the multi-grid method, which is faster than the conventional computation method (Murakami et al. 2023). We ran this solver on 40 NVIDIA A100 PCIe GPUs to achieve computation times as short as 2.9 h (10,500 s) per case of viscoelastic crustal deformation analysis.

Stress calculation on source faults

After obtaining the stress field σ_ij (x,t), we evaluated the normal and shear components of the stress on a fault plane as σ_ij (x,t) n_jn_i and σ_ij (x,t) n_jw_i, respectively, where n_i, and w_i denote the normal and slip vectors of the fault plane, respectively. As a metric of the preference of stress loaded on a fault toward rupture, the ΔCFS is widely used (Reasenberg and Simpson 1992; Stein et al. 1992). The ΔCFS is defined by the difference in shear stress and fault strength. The fault strength can be expressed by the normal compressive stress multiplied by the effective frictional coefficient μ′. Thus,

$$\Delta {\text{CFS}}\left( t \right) \, = \sigma_{ij} \left( {{\bf{x}},t} \right)n_j n_i + \mu^{\prime} \, \sigma_{ij} \left( {{\bf{x}},t} \right)n_j w_i .$$

(5)

Additional file 1: Table S1 shows the list of faults on which we evaluated the ΔCFS. We chose faults located west of 138°E in the Japanese islands from the fault list published by HERP (2009b). These faults are identified by geomorphological and geological techniques. However, HERP’s list does not cover all the source faults of the recent and historical damaging earthquakes. Thus, we added the source faults for the damaging earthquakes of interest, which are determined seismologically: the 1945 M6.8 Mikawa earthquake (Kakehi and Iwata 1992), the 2000 M7.3 west Tottori earthquake (Sagiya et al. 2002), and the 2004 M7.4 and M7.1 Kii Peninsula earthquakes (Suito 2017). We also add the primary nodal planes of the 1997 M6.6 Yamaguchi-Shimane earthquake, the 2001 M6.8 Geiyo earthquake, and the 2016 M6.6 central Tottori earthquake from the F-net catalog.

Results

Figure 4 shows the change over time of mean normal stress σ_mean = (σ₁₁ + σ ₂₂ + σ₃₃)/3 and von Mises stress σ_Mises at a depth of 10 km due to the 1944–1946 megathrust rupture sequence for Case 1 (uniform viscosity). For each stress component, we show the coseismic change at T₁₉₄₄ (Fig. 4a, b), viscoelastic change between T₁₉₄₄ and T₁₉₄₆ (Fig. 4c, d), coseismic change at T₁₉₄₆ (Fig. 4e, f), viscoelastic change between T₁₉₄₆ and T_End = T₁₉₄₄ + 4 (Fig. 4g, h), and total change in 4 years (Fig. 4i, j). Note that σ_mean and σ_Mises are related to the fault strength and the shear strain energy around the fault, respectively (Terakawa et al. 2020; Matsu’ura and Terakawa 2021).

Fig. 4

Change over time of the stress field at a depth of 10 km due to the 1944–1946 rupture for Case 1 (uniform viscosity). The left column shows the mean normal stress σ_mean, and the right column shows von Mises stress σ_Mises. a, b 1944 coseismic change, c, d 1944–1946 viscoelastic change, e, f 1946 coseismic change, g, h postseismic change after 1946, and i, j total stress change in 4 years

As shown in Fig. 4, the total stress change in 4 years is as large as ~ 0.1 MPa inland, while the viscoelastic stress change is 0.01 MPa at most near the source region. Thus, the inland stress is dominated by the 1944 and 1946 coseismic changes. Each rupture brings a ~ 0.1 MPa increase in σ_mean or decrease in fault strength, and also a change in σ_Mises, which significantly affects regional seismic activity.

To show the effect of viscoelastic relaxation on the stress change, we show viscoelastic stress change between T₁₉₄₄ and T₁₉₄₆ for Cases 1–5 (Fig. 5). Though the stress pattern due to the viscoelastic relaxation varies depending on the viscosity distribution, the magnitude of the stress change is too small, on the order of 0.001 MPa, to influence the overall stress change.

Fig. 5

1944–1946 viscoelastic change in σ_mean for five material cases

The viscoelastic effect is more visible in the stress change in the slab. Figure 6 shows the nearly 2-year change in σ_mean and σ_Mises under section BB′ from the 1946 rupture for Case 1. Additional file 1: Fig. S2 also shows the 4-year change in σ_mean and σ_Mises under section AA′ from the 1944 rupture for Case 1. Note that the effect on section AA′ from the 1946 rupture and that on section BB′ from the 1944 rupture are negligible. During each earthquake, the stress is distributed widely under its source region. As time passes, stress in the elastic slab is redistributed as viscoelastic relaxation progresses. As a result, the stress jump across the elastic–viscoelastic boundary, which is located 30 km below the plate boundary, becomes clearer. In contrast, the stress in the overriding plate shows almost no change, as we observe in Fig. 4.

Fig. 6

Change over time of the stress field due to the 1946 rupture under section BB′ for Case 1 (uniform viscosity). The left column shows the mean normal stress σ_mean, and the right column shows von Mises stress σ_Mises. a, b 1946 coseismic change, c, d viscoelastic change after 1946, and e, f total stress change in 2 years

Figure 7 and Additional file 1: Fig. S3 show the total stress change for Cases 1–5 in 2 and 4 years, respectively. In contrast to the inland stress fields (Fig. 4), we can see that the stress redistribution in the elastic slab is strongly influenced by the viscosity structure, particularly the existence of the LAB layer. In particular, the Case 4 result, which assumes an extremely low-viscosity LAB layer, shows complete decoupling of the stresses in the slab and the oceanic mantle in just 2 years. These results show the necessity of having knowledge of the viscosity structure for the seismogenesis in the slab in the postseismic stage.

Fig. 7

Total stress change under section BB′ for five material cases. The left column shows normal stress σ_mean, and the right column shows von Mises stress σ_Mises

Figure 8 shows the distribution of the ΔCFS on a map of the source faults. Here, the effective frictional coefficient μ′ is assumed to be 0.4 (King et al. 1994) to minimize uncertainty. Overall, the 1944 rupture makes positive ΔCFS on the strike-slip faults east of 135°E. In contrast, the 1946 rupture creates positive ΔCFS on the MTL and faults north of the MTL. As a result, the faults around 135°E, or the faults between the influential areas of the 1944 and 1946 ruptures, experience an inverse ΔCFS sign, from negative to positive. On Kyushu, the ΔCFS is generally negative through 4 years.

Fig. 8

ΔCFS values on source faults under SW Japan due to the 1944–1946 rupture sequence for Case 1 (uniform viscosity). The effective friction coefficient μ′ is set as 0.4. a ΔCFS at the 1944 rupture, b ΔCFS at the 1946 rupture, and c Total ΔCFS in 4 years. Contour of σ_Mises = 0.1 MPa is shown as influential areas of the 1944 and 1946 ruptures

To observe the 4-year stress change due to the 1944–1946 rupture for each fault in detail, we chose 11 faults from Fig. 1b that hosted M ≥ 6.8 earthquakes in the previous (1854–1944) and the current (1944-present) megathrust cycles. We classified these earthquakes by their location and time with respect to the 1854 or 1944 megathrust earthquakes as follows: (a) inland earthquakes after megathrust rupture (postseismic period), (b) inland earthquakes before megathrust rupture (interseismic period), and (c) intra-slab earthquakes. Here, we define the boundary of postseismic and interseismic periods to be 40 years after megathrust earthquakes so that the 1891 M8 Nobi earthquake, which occurred 37 years after the 1854 rupture, is included in Group a). For each group, we plot the 4-year change of stress for the end-member Cases 1 (uniform viscosity) and 4 (continental-oceanic viscosity contrast with an extremely low viscosity in LAB layer) (Figs. 9, 10, 11). Results for Cases 2, 3, and 5 are shown in Additional file 1: Figs. S4, S5, and S6, respectively. In each panel, we plot the following stress components: normal, shear, and ΔCFS for μ′ = 0.2, 0.4, and 0.6.

Fig. 9

Change over time of stress components for source faults for Group a) earthquakes (Table 2). In each panel, gray, red, green, blue, and purple lines denote normal stress; shear stress; and ΔCFS values for μ′ = 0.2, 0.4, and 0.6, respectively. The vertical line indicates the time of the 1946 rupture. The left column shows Case 1 (uniform viscosity), and the right column shows Case 4 (continental-oceanic viscosity contrast with an extremely low viscosity in LAB layer). a, b 1858 Hietsu earthquake, c, d 1891 Nobi earthquake, e, f 1945 Mikawa earthquake, and g, h 1948 Fukui earthquake

Table 2 List of damaging earthquakes for the source fault of which the 4-year stress change is shown in Figs. 9, 10, 11 and Additional file 1: Figs. S4–S6

Fig. 10

Same as Fig. 9 but for Group b) earthquakes. a, b 1927 North Tango earthquake, c, d 1995 Kobe earthquake, e, f 2005 Fukuoka earthquake, and g, h 2016 Kumamoto earthquake

Fig. 11

Same as Fig. 9 but for Group c) earthquakes. a, b 2001 Geiyo earthquake, c, d 2004 M7.4 Kii Peninsula earthquake, and e, f 2004 M7.1 Kii Peninsula earthquake

In Fig. 9 and Additional file 1: Fig. S4, we show the results for the Group a) faults, including the source faults of the 1858 M7.1 Hietsu and the 1891 M8.0 Nobi earthquakes, which occurred after the 1854 megathrust rupture, and those of the 1945 M6.8 Mikawa and the 1948 M7.1 Fukui earthquakes, which occurred after the 1944–1946 rupture. For these cases, the ΔCFS takes positive values, which is favorable for earthquake occurrence. As for the 1945 Mikawa earthquake, the ΔCFS amounts to several tens of kilopascals, which is consistent with the earthquake occurrence 1 month after the 1944 rupture. For the 1948 Fukui earthquake, the total ΔCFS was only several kilopascals. This value is very small compared to the stress drop value of ordinary earthquakes (on the order of megapascals). This stress change, however, corresponds to the steady stress accumulation over several to tens of years, considering the steady stressing rate in this region is nearly or less than 1 kPa (Hori and Oike 1999; Shikakura et al. 2014). On the source faults of the 1858 Hietsu and 1891 Nobi earthquakes, the ΔCFS takes positive values of several to 30 kilopascals. These example results suggest that the source faults of the postseismic inland earthquakes are similarly loaded in different megathrust cycles.

In Fig. 10 and Additional file 1: Fig. S5, we show the results for the Group b) faults, including the source faults of the 1927 M7.3 North Tango, 1995 M7.3 Kobe, 2005 M7.0 Fukuoka, and 2016 M7.2 Kumamoto earthquakes. For the source faults of the 2005 and 2016 earthquakes, the total ΔCFS is negative as other faults on Kyushu. Because stress on these faults decreases after the megathrust rupture, it takes time to recover the stress level prior to the megathrust rupture and to reach the fault strength. On Kyushu, M ≥ 6.5 inland earthquakes occurred in the interseismic period both in the previous and current megathrust cycles (Fig. 1c, d). In contrast to the above two earthquakes, the source fault of the 1995 Kobe earthquake shows a positive ΔCFS of tens of kilopascals and that of the 1927 earthquake show a positive ΔCFS of 1–2 kPa for all cases.

In Fig. 11 and Additional file 1: Fig. S6, we show the result for intra-slab earthquakes: the 2001 M6.8 Geiyo earthquake and the 2004 M7.4 and M7.1 Kii Peninsula earthquakes. The two earthquakes off the Kii Peninsula show a stress change as large as on the order of 100 kPa due to their proximity to the megathrust rupture, and the total ΔCFS become negative. In contrast, the 2001 Geiyo earthquake, which occurred at a depth of 50 km, shows only a 10 kPa change with a positive ΔCFS in 4 years. For these earthquakes, the contribution of the viscoelastic relaxation is also important. For the 2001 earthquake, the viscoelastic stress change is comparable to the coseismic change in Case 4. We can also observe a 10–100 kPa change in Case 4 for the Kii earthquakes.

Discussion

Framework of the postseimsic effect of the 1944–1946 Nankai trough ruptures

The occurrence of intra-plate earthquakes concentrates in a specific period and region within the cycles of large interplate earthquakes (Utsu 1974; Shimazaki 1976; Seno 1979; Ellsworth et al. 1981; Mogi 1981; Hori and Oike 1996; Bakun 1999). This is because the stress loading patterns in the plate interior due to locking and rupture on the plate interface are different. Intra-plate deformation due to the seismic cycle of large interplate earthquake was formulated by Savage (1983) and Matsu’ura and Sato (1989). According to them, interplate slip due to the seismic cycle can be represented by the superposition of steady uniform slip, an increase in the slip deficit due to locking, and repeating ruptures. The effect of slip deficit and repeating ruptures cancel in the long term. Only the effect of a steady uniform slip is permanently accumulated in the plate interior. However, it is much smaller in the short term than the effect of slip deficit and repeating ruptures. Of the contribution of all repeating ruptures, the effect of the latest rupture dominates, because an effect of a single megathrust rupture usually decays in the time scale of decades to hundred years. As slip deficit has been conventionally referred to as “back slip,” its effects are generally the inverse of rupture effects. Thus, the faults on which stress is accumulated or decreased after a megathrust rupture experience the lack of or accumulation of stress during the locking stage for the next megathrust rupture (Fig. 12a, b).

Fig. 12

a, b Schematic diagrams of two patterns of stress perturbation on an inland fault due to a Nankai megathrust cycle. Solid part of the curve denotes coseismic-postseismic stress, which is dominated by elastic and viscoelastic response of repeating ruptures. Broken part of the curve denotes interseismic stress, which is dominated by locking. Note that contribution from long-term steady sliding on the plate boundary is neglected. Shaded period denotes possible period of damaging earthquake occurrence. c, d Schematic diagrams of stress evolution in an inland earthquake cycle corresponding to stress perturbation pattern from the Nankai megathrust (a, b)

In the case of SW Japan, the inland region is not only affected by the seismic cycle of the Nankai trough, but also the steady east–west contraction due to the interaction between the PHS, Eurasia, and Pacific plates (Takahashi 2006, 2017; Fukahata 2019). Particularly in the inner zone, which is far from the Nankai trough, we need to take this effect into account. The effect of interaction between the intra-plate earthquakes is considered secondary (Hori and Oike 1999). The 1891 M8.0 Nobi earthquake exceptionally had a significant effect. The order of 100 kilopascals were loaded on the 1945 Mikawa and 1948 Fukui earthquakes (Pollitz and Sacks 1995; Shikakura et al. 2014). In the time period of our study, however, its viscoelastic effect is much smaller than the 1944–1946 sequence. The above-stated issues constrain the framework in discussing the relationship of stress accumulation models for SW Japan and intra-plate earthquake occurrence (Shikakura et al. 2014).

Figure 12 illustrates schematic diagrams of stress loading on an inland fault. Figure 12a, b show the close-up view of stress perturbation due to a Nankai megathrust cycle for the cases of the positive and negative coseismic loading, respectively. We show the dominant contributions from repeating ruptures and locking, neglecting the contribution from long-term steady sliding on the plate boundary. For an inland fault, the stress perturbation from a single megathrust cycle is not enough to reach the fault strength, which is typically of the order of ~ MPa. Superposing the long-term steady loading term, fault stress is gradually accumulated through multiple megathrust cycles to reach the fault strength, at which an inland earthquake occurs (Fig. 12c, d). In the last megathrust cycle before an inland earthquake, we can see that the stress is likely to reach the fault strength in a specific period; the postseismic period for the positive coseismic loading and the interseismic period for the negative coseismic loading due to a megathrust earthquake. This feature gives us a rough constraint on the period of possible inland earthquake occurrence on a fault by calculating coseismic stress change due to a megathrust earthquake even if we do not have information on the absolute stress level (Fig. 12a, b). Our stress change results in Figs. 9 and 10 on the source fault of the past damaging earthquakes are generally consistent with this view.

Focusing on the postseismic effect of megathrust earthquakes under the Nankai trough on the damaging earthquakes, we computed the stress field in SW Japan and stress loaded on source faults using the FE model which incorporates actual topography and elastic structure (JIVSM) around Japan. Studies on viscosity structures under plate subduction zones have rapidly accumulated in this decade by using the data from dense geodetic observation networks with improved analysis methods for postseismic deformation. However, we are still in the beginning stage of elucidating the deformation and viscosity structure. Based on Sun et al. (2014) and Suito (2017), we consider the viscosity of three regions: mantle wedge, oceanic mantle, and the LAB layer for the Nankai trough–SW Japan arc system. However, Freed et al. (2017) showed the depth (or temperature) dependence of viscosity in the mantle wedge and the oceanic mantle. We assumed 45 km for the thickness of the elastic continental lithosphere, but the proposed values vary from 30 to 60 km (e.g., Freed et al. 2017; Noda et al. 2018a,b). Also, the slab thickness assumed in this study might be underestimated considering the results of 60 km proposed by the seismic analysis of Kumar and Kawakatsu (2011) or geodetic analysis of Noda et al. (2018a,b). These effects should be examined in the future.

Viscoelastic relaxation

This study shows that viscoelastic relaxation from the megathrust events influences the intra-slab stress more than the inland stress in 4 years. This might be a natural result because the slab is closer to the source region and the LAB layer. In northeast Japan, M7 intra-slab earthquakes occurred in 2021 and 2022 in the postseismic stage of the 2011 M9 Tohoku earthquake, resulting in damage to the coastal area due to the strong ground motion. Zhang et al. (2023) showed that the 2021 intra-slab earthquake is likely promoted by the visoelastic relaxation with the ΔCFS increase of 0.3–0.7 MPa. To explain the mechanism of stress accumulation to intra-slab earthquakes, steady bending-unbending of the slab has been proposed (Christensen and Ruff 1988), but this effect should be unified with the effect of the seismic cycle (Matsu’ura and Sato 1989; Sato and Matsu’ura 1993). This study shows that intra-slab stress is particularly sensitive to the low viscosity of the LAB layer. Thus, we need to intensively study the viscosity of the LAB layer by updating past modeling studies with more accumulated geodetic data (Tomita et al. 2017; Yokota et al. 2018; Honsho et al. 2019).

As for the inland earthquakes, the viscoelastic effect seems negligible at least for the first 4 years because the coseismic effects of the 1944 and 1946 ruptures are dominant (Figs. 4, 5). For the inland region, the stress field would not be much different if we used the perfectly elastic model. Looking in detail at the temporal change in stress, however, viscoelastic relaxation might control the occurrence of some earthquakes depending on the choice of the viscosity structure. In the case of the source faults of the 1858 Hietsu and 1891 Nobi earthquakes, for example, the Case 4 result shows a decrease of the ΔCFS on the order of kilopascals over time after the megathrust rupture (Fig. 9b, d). If this trend continued, we cannot explain the 1891 earthquake occurrence 37 years after the 1854 megathrust rupture, provided that the viscoelastic trend of the 1854 rupture is similar to that of the 1944–1946 rupture and that the contribution from the steady east–west compression in the SW Japan arc are negligible. In contrast, the ΔCFS does not decrease for Cases 1, 2, 3, and 5. For these cases, the calculated stress does not contradict the earthquake occurrence. This example could be the constraint that viscosity in the LAB layer is not as low as the Case 4 value of 2.5 × 10¹⁷ Pa s. For more clear limitation of the viscosity in the LAB layer, we need more reliable seismicity data after the 1944–1946 rupture sequence particularly in the offshore.

Contributions from other processes

Postseismic deformation and stress are generally contributed by afterslip around the mainshock as well as viscoelastic relaxation in the asthenosphere. Because afterslip is usually localized in the deeper part of the mainshock, however, its contribution to the stress field is limited compared with that of viscoelastic relaxation (Becker et al. 2018). As for the 1944–1946 earthquake sequence, Sherrill and Johnson (2021) estimated the afterslip distribution separately from the viscoelastic relaxation. According to their results, the annual rate of afterslip is only 10 cm/year with the decay time of 15–30 years, which is on the order of 1/100 of the mainshock, even immediately after the mainshock. As our computation shows, the inland stress field is dominated by the coseismic effect of the 1944 and 1946 earthquakes. Thus, the effect of afterslip might be negligible.

We express the slip distribution of the 1944 and 1946 earthquakes by the superposition of five basic slip regions, referring to Sherrill and Johnson (2021) for each slip value. As we stated, estimates of the 1944–1946 slip show some differences depending on various previous studies. Figure 13 shows the ΔCFS patterns obtained by assuming different values for five basic slip regions of a = (3, 0, 8, 0, 8) and a = (2, 5, 0, 0, 12) by referring to Ichinose et al. (2003) for the 1944 rupture and Sagiya and Thatcher (1999) for the 1946 rupture, respectively, for Case 1. Though the absolute values of the ΔCFS differ on each fault, the overall characteristics shown in Fig. 8 hold; the positive ΔCFS is distributed in the eastern region due to the 1944 rupture and the positive ΔCFS occurs around the western MTL due to the 1946 rupture. In the middle of these regions near 135°E, we can see some difference due to the slip distribution. Particularly for the 1995 Kobe earthquake, the ΔCFS become positive after the 1946 rupture in Figs. 8 and 13a but negative in Fig. 13b, which references Sagiya and Thatcher (1999). Ogata (2004) investigated the relation between the activity of moderate earthquakes of M > 4 and the elastic stress change due to the 1944 and 1946 megathrust earthquakes, concluding that the calculated positive ΔCFS is consistent with these seismic activities (Figs. 5F, 5I, 5J, and 5K in Ogata 2004, corresponding to our result). However, the seismic activity in Ogata includes the aftershocks of local M7 earthquakes. Comparing the 20-year activity before the 1944 rupture and activity after the 1946 rupture, we cannot conclude that the seismic activity increased, particularly for the 1927 and 1995 earthquakes (Figs. 5K and 5F in Ogata 2004).

Fig. 13

Total ΔCFS in 4 years on inland faults under SW Japan due to the 1944–1946 rupture sequence and change over time of stress components for source faults for the 1995 Kobe earthquake for different slip settings from Fig. 9 for Case 1 (uniform viscosity). a, b a = (3, 0, 8, 0, 8) by referring to Ichinose et al. (2003) for the 1944 rupture, and c, d a = (2, 5, 0, 0, 12) by referring to Sagiya and Thatcher (1999) for the 1946 rupture

If the Fig. 13b result is correct, the 1995 earthquake was caused by a stress loading process similar to that of the 2005 Fukuoka and 2016 Kumamoto earthquakes, where the stress decrease after the megathrust rupture recovered and further accumulated to the level of the fault strength due to the steady stress loading (Fig. 10e–h). This result implies a preference for the 1995 earthquake occurrence explanation of Sagiya and Thatcher (1999), where slip concentrates under the Tosa Bay in the westernmost slip region, over that of Sherrill and Johnson (2021), where slip widely distributes along the trench axis. Sherrill and Johnson (2021) adopted a two-dimensional plate boundary model, but the actual plate boundary largely bends around 135°E. This difference is likely to strongly affect the result of Sherrill and Johnson (2021). As shown in the above example, stress computation compared to the seismicity might constrain the slip distribution for the different slip models.

This study compared the calculated stress with the M ≥ 6.5 earthquake occurrence, including those that occurred before modern observations. Use of large earthquakes may lack statistical robustness because of their scarcity. The analysis of Ogata (2004) included moderate earthquakes of M ≥ 4. However, small earthquakes have the nature of clustering in time and space as the aftershocks of larger events. It is preferable to use declustered seismicity for the comparison. In addition, for smaller earthquakes, it is difficult to determine the fault plane with accuracy and to directly apply the ΔCFS to them. Recently, Terakawa et al (2020) proposed a new metric called the change in the energetic-based failure stress, which is a generalized metric of the ΔCFS. This new metric is promising for elucidating the activity of smaller earthquakes. A unified understanding of historical damaging earthquakes and modern smaller earthquakes would be an issue for future work.

Another approach in the future may be to physically reveal how modeled stress from the megathrust rupture nucleates large earthquakes on inland and intra-slab faults. Parsons (2002) and Gomberg et al. (2005) explained the Omori law of aftershock decay by implementing rate- and state-dependent friction (Dieterich 1994). Cascading nature of the dynamic rupture starting from the nucleation phase (e.g., Ellsworth and Beroza 1995) is also important. This effect may work differently on the seismicity depending on the style of the stress loading (stepwise or gradual loading). Combining this effect with our viscoelastic stress loading model may reveal the physical nature of the inland earthquake occurrence due to the 1944–1946 sequence.

Conclusions

Using a high-fidelity FE model incorporating actual topography and a plausible viscoelastic structure, we estimated the 4-year stress change in SW Japan due to the 1944 Tonankai and 1946 Nankai earthquakes. Our calculated stress field shows the dominance of coseismic change in 1944 and 1946 and also the strong effect of viscoelastic change in the slab, indicating the importance of quantifying viscosity structures, particularly the LAB layer.

The ΔCFS values are positive on strike-slip faults east of 135°E due to the 1944 rupture and on the MTL and faults north of the MTL due to the 1946 rupture. As a result, in the middle of the influential region of the 1944 and 1946 ruptures, we can observe inversion of the sign of the ΔCFS from negative to positive. The faults in Kyushu maintained a negative ΔCFS in our calculation period. These trends can explain the occurrence of damaging inland earthquakes such as the 1945 Mikawa and 1948 Fukui earthquakes. The negative ΔCFS on the source faults of interseismic earthquakes such as the 2016 M7.2 Kumamoto earthquakes shows the delay of stress accumulation. These trends do not depend much on the viscosity structure in general. However, the postseismic ΔCFS on the source fault of the 1891 Nobi earthquake differs slightly depending on the viscosity in the LAB layer. Another exception is the 1995 Kobe earthquake, the analysis of which was sensitive to the slip distribution of the 1944–1946 sequence rather than the viscosity structure. These results might be useful for establishing constraints on the viscosity structure or the slip distribution. On the faults of the intra-slab earthquakes, the ΔCFS depends significantly on the viscosity of the LAB layer.

In the future, the slip distribution used in this study should be updated by the inversion using the FE-model-based Green’s functions for performing a more realistic stress computation. This would be useful for revealing the mechanism of damaging inland earthquakes and the viscosity structure by comparing the result with the historical damaging earthquakes and modern smaller earthquakes in a unified way. For this purpose, the geometry of the source faults and the rheological properties around them should be also investigated with more accuracy.