Coulomb stress change on inland faults during megathrust earthquake cycle in southwest Japan

In the subduction zone, megathrust earthquakes may modulate the shallow crustal seismicity in the overriding plate. Historical documents indicate the frequent occurrence of large shallow crustal earthquakes in the overriding continental plate 50 years before and 10 years after the megathrust earthquakes along the Nankai trough in southwest Japan. In this study, we model megathrust earthquake cycles in a simple oblique subduction zone considering the viscoelasticity, and calculate the temporal evolution of the Coulomb failure stress changes (ΔCFS\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Delta {\text{CFS}}$$\end{document}) on the crustal faults in the overriding plate. Further, we examine the variation of ΔCFS\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Delta {\text{CFS}}$$\end{document} depending on the location and fault type, and the active period of crustal earthquakes in which ΔCFS\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Delta {\text{CFS}}$$\end{document} exceeds the previous maximum. Our viscoelastic model suggests that the dependency of the active period on the distance from the megathrust fault is less when the intrinsic loading rate of the inland fault is low. Moreover, it suggests that the viscoelastic stress evolution on faults with negative coseismic ΔCFS\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Delta {\text{CFS}}$$\end{document} renders the active period longer or shorter than those in a pure elastic medium. The temporal evolution of ΔCFS\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Delta {\text{CFS}}$$\end{document} on most major active faults in southwest Japan can be categorized into two groups with the following different characteristics: one is that ΔCFS\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Delta {\text{CFS}}$$\end{document} is positive coseismically and peaks 10 years after a megathrust earthquake. The other is that ΔCFS\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Delta {\text{CFS}}$$\end{document} is negative coseismically, and does not recover to the preseismic one for more than 50 years after a megathrust earthquake. This can explain the temporal sequence of the historical earthquakes in southwest Japan. Our model which includes viscoelastic relaxation successfully expresses the activation of shallow crustal earthquakes in the overriding continental plate not only before the megathrust earthquake, but also after. If the apparent frictional coefficient is less than ~ 0.1, the coseismic ΔCFS\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Delta {\text{CFS}}$$\end{document} on the source faults of the 1943 Mj7.3 Tottori earthquake, 1596 M7.0 Keicho Iyo earthquake, and 1596 M7.5 Keicho Fushimi earthquake that occurred within 10 years before the megathrust earthquake along the Nankai trough is negative. Therefore, to explain the occurrence of these historical earthquakes, our model suggests that the apparent frictional coefficient must be less than ~ 0.1.


Introduction
Megathrust earthquake cycles in the subduction zone can modulate the shallow crustal seismic activity in the overriding continental plate through quasi-static stress changes. Megathrust earthquakes occur at a 100-200 year intervals with the subduction of the Philippine Sea Plate in southwest Japan (Earthquake Research Committee 2001). Previous studies (e.g., Ustu 1974;Hori and Oike 1996) indicate the frequent occurrence of shallow crustal earthquakes in the overriding continental plate (hereafter, inland earthquakes) 50 years before and 10 years after the great megathrust earthquakes along the Nankai trough (Fig. 1). Southwest Japan is one of the best regions to examine the temporal relation between inland earthquakes and megathrust earthquakes due to the availability of the longest, and most continuous and complete set of historical records of large earthquakes in the world and detailed datasets for the distribution of inland Open Access *Correspondence: mitogawa.tsukasa.53w@st.kyoto-u.ac.jp 1 Graduate School of Science, Kyoto University, Uji 611-0011, Japan Full list of author information is available at the end of the article Page 2 of 10 Mitogawa and Nishimura Earth, Planets and Space (2020) Utsu (1990Utsu ( , 2002Utsu ( , 2004, whereas that from 1922 to 2018 is based on the Japan Meteorological Agency (JMA). The epicenters are plotted north of the black solid line in the range of a. The light and dark red circles indicate the epicenters of crustal earthquakes that occurred within 50 years before and 10 years after a megathrust earthquake. The light and dark blue circles indicate the epicenters in the other periods. The light and dark colored circles indicate the epicenters during 697-1921 and 1922-2018, respectively active faults (e.g., Utsu 1990Utsu , 2002Utsu , 2004Earthquake Research Committee 2018). Many previous studies have quantitatively evaluated the effect of megathrust earthquakes on the occurrence of inland earthquakes (e.g., Pollitz and Sacks 1997;Hori and Oike 1999;Shikakura et al. 2014). An earthquake generally occurs when shear stress on the fault exceeds its strength. Although some studies including the fault valve model (Sibson 1990) consider temporal evolution of the fault strength, many models focus on stress evolution on the fault to examine when the earthquake occurs (e.g., Pollitz and Sacks 1997;Shikakura et al. 2014). Coulomb failure stress ( CFS ) and its change ( CFS ) are often used for evaluating stress on the fault and have been applied for evaluating the modulation of the inland seismicity in the overriding plate by a megathrust fault. Hori and Oike (1999) applied the elastic response function for calculating CFS on the crustal fault by megathrust earthquakes, and numerous studies applied the viscoelastic response function (e.g., Pollitz and Sacks 1997;Shikakura et al. 2014). Pollitz and Sacks (1997) succeeded in explaining the occurrence of the 1995 M j 7.3 Kobe earthquake due to delayed CFS increase after the megathrust earthquakes in 1944 and 1946, using the viscoelastic response function. Shikakura et al. (2014) succeeded in explaining most of the inland earthquakes in the Kinki region, incorporating the interaction between inland earthquakes, using the viscoelastic function.
Although previous studies developed specific models to focus on reproducing these historical events, the general tendency of the inland fault response to interplate earthquakes and megathrust fault locking has not been investigated. In view of the above, we model the interplate earthquake cycles for a simple oblique subduction zone and calculate the temporal evolution of CFS on inland faults in this study. We examine the variation of CFS depending on the fault type, intrinsic loading rates, and distance from the megathrust. Furthermore, we apply the calculated CFS to the inland faults of southwest Japan, and discuss the general influence tendency by a megathrust earthquake cycle.

Method
Coulomb failure function ( CFS ) is defined as where τ s , σ n , and µ ′ are the shear-stress, normal-stress (positive in extension) with the fault geometry (strike dip and rake), and the apparent friction coefficient, respectively. Therefore, the stress change of CFS ( CFS ) is defined as where �τ s and �σ n are the changes in the shear-stress, normal-stress with the fault geometry, respectively. Following many previous studies using CFS (e.g., Pollitz and Sacks 1997;Hori and Oike 1999;Shikakura et al. 2014), the apparent friction coefficient is assumed to be constant. The apparent friction coefficient is intended to include the effects of pore fluids as well as the material properties of the fault zone (Harris 1998). At 2nd line of Eq.
(2), we divided �τ s and �σ n into external perturbation components from a megathrust fault (i.e., �τ external s and �σ external n ) and secular intrinsic loading ones on the inland fault (i.e., �τ intrinsic s and �σ intrinsic n ), respectively. The former process comprises interseismic locking and coseismic slip in the megathrust fault, whereas the latter can arise from background tectonic loading including a regional stress increase. Although we support aseismic sliding on the downward extensions of the seismogenic part of the inland fault (Iio and Kobayashi 2002) for the source of intrinsic loading, we do not specify the source and simply assume a constant shear loading rate with time for all the faults. In other words, we assigned a certain value of �τ intrinsic s per year and �σ intrinsic n = 0 for the intrinsic loading except for the case of the regional eastwest compression in "Discussion" section.
External CFS perturbations on inland faults are calculated using the viscoelastic response function for explaining the delayed stress change response after a megathrust earthquake, in particular. The stress change due to a megathrust fault was computed using numerical code by Fukahata and Matsu'ura (2006), assuming a medium composed of an elastic layer overlying a Maxwell viscoelastic half-space as shown in Additional file 1: Fig. S1. The effect of gravity was incorporated as a buoyant force caused by vertical displacement on the earth's surface. The thickness of the elastic layer was 30 km. The density, bulk modulus, and shear modulus were 3.00 × 10 3 kg/ m 3 , 58.4 GPa, and 40.0 GPa, respectively, in the elastic lithosphere, whereas the density, bulk modulus, shear modulus, and viscosity were 3.40 × 10 3 kg/m 3 , 130 GPa, 60.0 GPa, and 1.00 × 10 19 Pa s in the viscoelastic substratum, respectively. The Maxwell relaxation time of the viscoelastic material was 5.28 years.
We focus on CFS depending on the distance from the megathrust fault and do not consider the along-strike variation, i.e., our model setting simulates the stress change across an infinitely long megathrust fault. In practice, we approximate it by a fault length of 2000 km, and calculate the stress change across the center of the megathrust fault. For simplicity, we calculate the stress change at a depth of 10 km for all the inland faults. The megathrust fault model has a relatively low dip angle and oblique subduction zone, approximating the subducting Philippine Sea Plate along the Nankai trough. The megathrust fault is represented by a single rectangular fault whose depth of the upper edge, width, strike, dip, and rake are 11 km, 108 km, 180°, 10°, and 117°, respectively. A coseismic slip of 6 m is provided every 100 years. The interplate locking in the interseismic period can be expressed as the superposition of the steady plate subduction over the entire plate interface with a back-slip in the locked portion (Savage 1983). Because the shape of the plate interface is flat, we neglect the effect of steady plate subduction. The back-slip rate is assumed to be 6 cm/year (= 6 m/100 year), and constant during the interseismic period such that all the stress accumulated during the interseismic period are released by the coseismic slip. The stress loading rate due to steady locking can be expressed by the complete relaxed response to an instantaneous back-slip in the viscoelastic medium (e.g., Matsu'ura and Sato 1989). We apply the response after 10,000 years of slip as the response due to steady locking. We repeated 18 megathrust earthquake cycles until the stress change caused by the megathrust fault became approximately cyclic after the initial effect. We present the stress evolution during the 19th and subsequent earthquake cycles in the next chapter. Further, we investigate the temporal evolution of CFS with varying the apparent friction coefficient and the intrinsic loading rate.

Coseismic stress change
We calculated the coseismic stress change on the inland fault at a distance of 150 km from the megathrust fault.
In the case of a strike-slip fault, the polarity of the coseismic-shear-stress change is highly dependent on the strike of the inland fault. The coseismic-shear-stress change on a reverse fault is negative with wider range of strikes and dips, whereas that on a normal fault is positive (Fig. 2a). It is notable that a larger value of µ ′ generally increases CFS because the coseismic-normal-stress change is positive with wider range of strikes and dips (Fig. 2b). Thus, the coseismic polarity of CFS on inland faults depends on the value of µ ′ . For example, the nor-

Temporal CFS evolution incorporating the viscoelasticity
Viscoelasticity incorporation causes stress-rate changes during the interseismic period. In addition, it results in different interseismic stress-change characteristics depending on the distance from the megathrust fault. For example, for the dextral and sinistral faults located 150 km and 300 km from the upper-edge of the megathrust fault, the used strike, dip, intrinsic loading rate ( �τ intrinsic s /year), and apparent friction coefficient are 120°, 90°, 1.0 kPa/year, and 0.1, respectively (Fig. 3). These distances from the megathrust upper-edge roughly correspond to the actual range of the major active faults in the studied area (Fig. 1). CFS increases immediately after the coseismic CFS drop caused by a megathrust earthquake in case of the dextral fault locates at a distance of 150 km (solid line in Fig. 3a). In contrast, the postseismic CFS decreases for a duration after the coseismic CFS drop at the distance of 300 km (dotted line in Fig. 3a).
If the coseismic CFS is positive (Fig. 3b), CFS exceeds the previous maximum when a megathrust earthquake occurs. When sinistral inland faults are located 150 km from the megathrust fault, their CFS exceeds the previous maximum only at the time of the earthquake (solid line in Fig. 3b); however, in case they are located 300 km from the megathrust fault, their CFS progressively exceeds the previous maximum for 33 years after the megathrust earthquake (dotted line in Fig. 3b). On the other hand, CFS increases during the interseismic period in the case of negative coseismic CFS ; it exceeds the previous maximum after a considerable period following a megathrust earthquake, and continues to exceed until the next megathrust earthquake. Examples of dextral faults at the distance are 150 and 300 km from the distance from the megathrust fault (Fig. 3a) indicate that CFS exceeds the previous maximum 30 years and 46 years before the megathrust earthquake, respectively. Although we know neither absolute value of CFS nor fault strength in terms of earthquake occurrence, it is certain that an earthquake cannot occur as long as CFS is below its previous maximum. Therefore, an earthquake can occur only during the period when CFS exceeds its previous maximum and we defined this period as 'the active period' of inland earthquakes (Fig. 3). The active period varies depending on the intrinsic loading rate of the inland fault. We investigate the active period variation with the distance and the intrinsic loading rate for four different fault types of inland earthquakes (Fig. 4). Here, negative and positive values of the active period represent activated inland seismicity before and after the megathrust earthquake, respectively. When the distance from the megathrust fault is small (150-200 km), the active period after the megathrust earthquake is short because CFS decreases immediately after the earthquake. When the distance from the megathrust fault is longer than 200 km, the active period after the megathrust earthquake is longer and depends on the intrinsic loading rate. On the other hand, the active period before the megathrust earthquake changes considerably depending on the intrinsic loading rate, even when the distance from the megathrust fault is less. It is interesting that the dependency of the active period on the distance from the fault is less when the intrinsic loading rate is low.

Impact of viscoelastic relaxation on the active period
If the viscoelasticity is incorporated, the variation of the interseismic stress rate of a fault with negative coseismic CFS differs depending on the distance from the megathrust fault (Fig. 3). In a pure elastic medium, because CFS increases at a constant rate during the interseismic period of a megathrust earthquake, the active period can be easily calculated. For example, the active period is 13 years and 30 years for the pure elastic and viscoelastic a b Fig. 3 Examples of the temporal evolution of CFS on an inland fault. The apparent coefficient ( µ ′ ) and the intrinsic loading rate ( �τ intrinsic s /year) in the inland fault are 0.1 and 0.5 kPa/year, respectively. a Case of negative coseismic CFS for dextral faults. The solid and dotted blue lines indicate CFS for faults located 150 km and 300 km from the upper-edge of the megathrust fault, respectively. The dark and light blue-shaded areas represent the active periods of these faults, respectively, and b case of positive coseismic CFS for the sinistral faults. The solid and dotted red lines indicate CFS for faults located 150 km and 300 km from the upper-edge of the megathrust fault, respectively. The dark and light red-shaded areas represent the active periods of these faults, respectively. However, the dark red-shaded areas show only during megathrust earthquakes Mitogawa and Nishimura Earth, Planets and Space (2020) 72:66 cases, respectively, when the intrinsic loading rate is 1.0 kPa/year for a vertical dextral fault whose distance from the upper-edge of the megathrust fault and strike are 150 km and 120°, respectively (Additional file 1: Fig.  S2a). The active period is 68 years and 46 years for the pure elastic and viscoelastic cases, respectively, with the same fault mechanism but at a distance of 300 km (Additional file 1: Fig. S2b). The viscoelastic effect on the active period of inland earthquakes is generally not straightforward, and renders the active period longer or shorter than that of a pure elastic medium depending on the intrinsic loading rate, focal mechanism, and distance from the megathrust fault.

Application to major active faults in southwest Japan
The obtained results mentioned in the previous chapter were applied to the major inland active faults in southwest Japan (Earthquake Research Committee 2018). The location, geometry, and rake of the inland active faults were referenced from the Earthquake Research Committee (2018). As our model does not consider the variations along the strike of the megathrust fault, we set the upper-edge of the megathrust fault with a strike of 241°, as denoted by the solid line in Fig. 5. We used two values of the apparent frictional coefficient ( µ ′ = 0.1 and 0.4 ), because it can affect the polarity of the coseismic CFS . For µ ′ = 0.4 , most of the reverse faults in the Median Tectonic Line (MTL) in the eastern Shikoku and Kinki districts exhibited a tendency to be activated before the megathrust earthquake. On the other hand, the dextral faults in the Kinki, northern Chugoku, and western Shikoku districts were activated after the earthquake; however, some of the dextral faults were activated before the megathrust earthquake for µ ′ = 0. Keicho Fushimi earthquake, respectively (Fig. 5). These large inland earthquakes occurred within 10 years before the 1605 M7.9 Keicho Nankai earthquake and the 1944 M j 7.9 Showa Tonankai earthquake, which are megathrust earthquakes along the Nankai trough. Therefore, we propose the possibility that the apparent frictional coefficient was less than ~ 0.1 to explain the occurrence of these historical earthquakes. Our inference on the small frictional coefficient is consistent with those of previous studies. For example, Shikakura et al. (2014) suggested that the apparent frictional coefficient of major inland active faults may be less than 0.15 because CFS during the interseismic period of inland earthquakes for certain faults is negative due to tectonic east-west compression, if the apparent frictional coefficient is 0.15 or more. Iio (1997) estimates that the friction coefficient associated with a fault in a seismogenic region is smaller than 0.2 from the numerous focal mechanism solutions in the Kinki region.
The friction coefficient of rocks has been shown to be about 0.65-0.8 for most rocks (Byerlee 1978). However, friction coefficients for many fault gouge minerals in wet condition can be dropped by 20-60% from the dry condition (Morrow et al. 2000). In southwest Japan, a b c d  Nishimura Earth, Planets and Space (2020) 72:66 dehydration from the subducted oceanic slab occurs due to various processes including diagenesis and denaturation. The dehydrate water from the slab may drop the apparent friction coefficient of the fault zone in the inland arc crust. Although we applied a uniform stress rate for all the inland faults in this study, the intrinsic loading rate actually differs for each inland fault. Here, we test a case of the uniform regional stress similar to the case of Shikakura et al. (2014). When we assume the uniaxial compressional stress rate of 2.9 kPa/year comparable to geological and seismological strain rates ( ∼ 3 × 10 −8 / year) in a direction of N100˚E, the intrinsic CFS loading rate on the inland faults is calculated to be between − 0.5 and 1.5 kPa/year when the apparent friction coefficient is 0.1 (Additional file 1: Fig. S3c). The spatial pattern of the active period (Additional file 1: Fig. S4) is similar to those in a case of a uniform loading rate of 1.0 kPa/year (Fig. 5a). On the other hand, a strain rate observed by geodetic measurements is about ~ 3-5 times larger than rates estimated by geological and seismological studies in Japan (Shen-Tu et al. 1995). When the intrinsic loading rate is 5 kPa/year, the simulated active period is long (Fig. 5b, d) and an inland earthquake can occur at any time during the megathrust earthquake cycle. In order to explain the modulation of historical earthquakes, the loading rate must be much smaller than that simply converted from the geodetic strain rate using linear elasticity.
Our assumption on the intrinsic loading, that is, a uniform rate for all faults and a uniform regional stress, is probably an unrealistic oversimplification. The uniform regional stress contradicts with large variations on longterm slip rates and recurrence intervals of earthquakes among inland faults with similar fault geometry in southwest Japan. GNSS observation revealed highly heterogeneous distribution of strain rates in southwest Japan (Sagiya et al. 2000;Nishimura et al. 2018) and it should be possible to determine the current stress loading rate of individual inland faults by incorporating the observed crustal movements such as GNSS data. Moreover, it is necessary to consider the interaction between inland earthquakes, particularly in areas where the active inland faults are densely distributed as in the Kinki district. These subjects are intended for future research.

Conclusions
We investigated the influence of megathrust earthquake cycles on inland faults in the oblique subduction zone using the temporal evolution of CFS considering the viscoelasticity. In order to focus on CFS depending on the distance from the megathrust fault, the megathrust fault length was set semi-infinite in the strike direction, in our model. CFS on inland faults are caused by megathrust earthquakes and locking, and the intrinsic loading on the inland fault.
Our model indicated that the polarity of the coseismic CFS is highly sensitive to the inland fault strike in the case of a strike-slip fault, but not in the case of a dip-slip fault. The polarity of the interseismic CFS is opposite to that of the coseismic CFS , without considering the viscoelastic effect. However, the polarity of CFS may change to that of the coseismic one during an interseismic period, if the viscoelasticity is incorporated. As a result, the viscoelastic stress evolution on faults with negative coseismic CFS can render the active period of inland earthquakes longer or shorter than those of a pure elastic medium. On applying the obtained results to active faults in southwest Japan, it was determined that certain active faults which are the source faults of historical earthquakes that occurred before the megathrust earthquake can exceed the maximum CFS before a megathrust earthquake only when the apparent friction coefficient is less than ~ 0.1. Therefore, it was suggested that the apparent friction coefficient of these inland faults is less than ~ 0.1. The low apparent friction coefficient in southwest Japan may indicate weak material properties of the fault zone.