Earthquake triggering model based on normal-stress-dependent Nagata law: application to the 2016 Mie offshore earthquake

We propose a normal-stress-dependent Nagata law. Nagata et al. (J Geophys Res 117:B02314, 2012) revised the rate- and state-dependent friction law by introducing the shear stress dependence. We further extended the Nagata law by incorporating the normal stress dependence obtained by Linker and Dieterich (J Geophys Res 97:4923–4940, 1992). We performed numerical simulations of earthquake triggering by assuming the extended Nagata law. In the case of repeated earthquakes, we applied dynamic Coulomb failure function (CFF) perturbation due to normal or shear stress changes. CFF perturbation increased the slip velocity after the cessation of perturbation, relative to that of the repeated events without triggering. This leads to dynamic earthquake triggering for certain perturbation amplitudes with time to instability of 0 to several tens of days. In addition, triggering potential of the static CFF jump (ΔCFFs) was investigated. Static stress perturbation has a higher triggering potential than dynamic stress perturbation for the same magnitude of CFF. The equivalent ΔCFFeq is evaluated for dynamic perturbation that results in a triggering potential approximately the same as in the case of static stress perturbation if ΔCFFs = ΔCFFeq. We calculated ΔCFFeq on the interface of the Philippine Sea plate for the Mie offshore earthquake, which occurred around the Nankai Trough on April 1, 2016, using OpenSWPC. The results shows that ΔCFFeq is large around the trough, where slow slip events followed the Mie earthquake, suggesting that a large ΔCFFeq may have triggered slow slip events.


Introduction
There have been many reports of earthquake triggering (e.g., Hill et al. 1993;Brodsky et al. 2000;Kilb et al. 2000;Miyazawa et al. 2008;van der Elst and Brodsky 2010;Richards-Dinger et al. 2010;Miyazawa 2011;Parsons and Velasco 2011;Lin 2012;Pollitz et al. 2012;van der Elst et al. 2013;Yukutake et al. 2013;Johnson and Bürgmann 2016;Miyazawa 2016;Uchida et al. 2016;Yoshida 2016). Besides the seismic event triggering, triggering of very low-frequency earthquakes, tremors, slow slip events (SSEs), and creep events has been observed (e.g., Allen et al. 1972;Obara 2009;Itaba and Ando 2011;To et al. 2015;Araki et al. 2017;Miyazawa 2019;Canitano et al. 2019;Katakami et al. 2020). A likely mechanism behind earthquake triggering has been based on the rate-and state-dependent friction (RSF) law (Dieterich 1979;Ruina 1983). Dieterich (1994) obtained the time to instability as a function of slip velocity based on a spring-block model, and found that an increase in the slip velocity due to a static stress change could trigger an earthquake. This model has been extended to include dynamic stress perturbation, which reduces the friction strength (expressed as a state variable), causing an increase in slip velocity (e.g., Yoshida 2018).
Numerical simulations of dynamic triggering have been extensively performed by Gomberg et al. (1997Gomberg et al. ( , 1998, Belardinelli et al. (2003), Savage and Marone (2008), and van der Elst and Savage (2015) using a spring-block model, and by Perfettini et al. (2003), and Kaneko and Lapusta (2008) using a one-dimensional (1D) fault in a two-dimensional (2D) medium in the RSF framework. Besides the seismic event simulation, Du et al. (2003), Perfettini and Ampuero (2008), and Wei et al. (2015) performed numerical simulations of triggering of SSEs, assuming velocity weakening region with a system stiffness smaller than but close to a critical value, or velocity strengthening region. Although the slip velocity acceleration after the start of the instability is different between the seismic event and SSEs, the growth of the rupture nucleation leading to the instability is basically the same, and hence, the same triggering mechanism is applicable. Yoshida (2018) performed numerical simulations of earthquake triggering by assuming a circular asperity that obeys a RSF law revised by Nagata et al. (2012). He also discussed delayed and instantaneous triggering due to changes in dynamic and/or static stress assuming a constant normal stress. In the present study, we will incorporate the normal stress dependence obtained through normal stress jump experiments by Linker and Dieterich (1992) into the Nagata law, and discuss Coulomb stress change (CFF) for this law in "Normal-stress-dependent Nagata law" section. With regard to dynamic stress perturbation, we will define ΔCFF eq so that the triggering potential of the dynamic perturbation is approximately the same as that of the static CFF change if ΔCFF eq is equal to the value of static CFF. In "Simulation of earthquake triggering" section, we perform a numerical simulation of earthquake triggering. In "CFF due to the 2016 Mie offshore earthquakes" section, we will calculate ΔCFF eq on the interface of the Philippine Sea plate due to the Mie offshore earthquake, which occurred around the Nankai Trough on April 1, 2016.

Normal-stress-dependent Nagata law
On the basis of the friction experiments for constant normal stress, Nagata et al. (2012) provided a revised RSF law, which is given as: where τ is the shear stress, σ is the effective normal stress, a and b are friction parameters, L is the characteristic slip, θ is a state variable, V* is an arbitrarily chosen reference velocity, θ* is a reference state at V = V*, and c is the stress-weakening effect coefficient. The slip law (Ruina 1983) and the aging law (Dieterich 1979) are used most frequently in numerical simulations. However, the former cannot reproduce the log t healing at zero slip velocity, and the latter cannot represent symmetric stress changes for sudden increases and decreases of the slip velocity. The Nagata law can reproduce both of these behaviors found in the laboratory experiments. As τ = μσ, Eq.
(2) may be written as: which apparently represents normal stress dependence. However, Eq.
(2) was derived from the experimental data for constant normal stress, so that the third term cµθ bσ dσ dt on the right-hand side of Eq. (3) was not experimentally confirmed. Linker and Dieterich (1992) conducted normal stress step experiments, keeping constant sip velocity V, and proposed the following evolution law that considers the effect of variable normal stress: where α is a constitutive parameter representing normal stress dependence. Incorporating the third term of Eq. (4) to Eq. (2), we propose the following normalstress-dependent Nagata law as: which includes laboratory-derived normal and shear stress dependences. Although Eq. (3) can be regarded as a specified version of Eq. (5) for α = cμ, since Nagata et al. (2012) derived Eq. (2) on the basis of a laboratory experiment under the constant normal stress, they did not check whether the normal stress change followed cµθ bσ dσ dt or not.
From Eq. (5) and the time derivative of Eq. (1), we obtain: , we obtain: Linker and Dieterich (1992) performed normal stress step tests where the normal stress was suddenly changed from σ 0 to σ 1 = σ 0 + Δσ keeping the constant sliding . Yoshida et al. Earth, Planets and Space (2020) 72:141 velocity. In such a case, the change in ln(V/V * ) due to the normal stress jump is written as: where Equation (9) defines ΔCFF for the extended Nagata law, which coincides with usual ΔCFF when α = c = 0. The lefthand side of Eq. (8) is zero because of the constant V, and the second term on the right-hand side is approximately zero because Δt takes a very small value and θ has a steadystate value of L/V immediately before the normal stress jump so that ΔCFF is also approximately zero. Therefore, the increase of the shear stress is written as: Linker and Dieterich (1992) defined δ L τ as the difference between the stress immediately after the stress jump and the following steady-state stress, as shown in Fig. 1a. For the case of the normal-stress-dependent Nagata law, Eq. (10) leads to Linker and Dieterich (1992) found that δ L μ = δ L τ/σ changes as follows: Comparing Eqs. (11) and (12), scaling relation for α between the normal stress-dependent aging law and the Nagata law is represented as: Figure 1a shows the simulation result of the normal stress step test from σ 0 = 5 MPa to σ 1 = 6 MPa using the 1D spring-block model for the aging law (4). The assumed parameters are a = 0.012, a-b = − 0.0015, L = 2.25 × 10 −6 m, α = 0.2, and V = 1 × 10 −6 m/s. Figure 1b shows the result when assuming the normal stress-dependent Nagata law (5). To reproduce the change in the same way as the aging law, the following scaling relationships were given by Bhattacharya and Rubin (2014) as: The friction parameters were converted using the scaling relationships (13) and (14) with c = 1. The shear stress for the Nagata law changes in a similar way as that for the aging law. Figure 1c shows the simulated relationship between δ L μ and ln(σ 1 /σ 0 ) for the aging law and (14) Fig. 1 Simulation of the normal stress step test. During stable sliding at a constant velocity, we applied static normal stress change σ 0 = 5 MPa to σ 1 = 6 MPa. a Shear stress change assuming the aging law. a = 0.012, a−b = − 0.0015, L = 2.25 × 10 −6 m, and α = 0.2. b Shear stress change assuming the normal-stress-dependent Nagata law. The friction parameters were converted using the scaling relations (14) and (15) with c = 1. c The relationship between δ L μ and ln(σ 1 /σ 0 ) for the aging law and the Nagata law the Nagata law. Indeed, the Nagata law with the scaling relationships reproduces frictional changes similar to those derived from the aging law. In these simulations, we assumed a very high stiffness of the slider to keep the sliding velocity constant at the time of the sudden normal stress change, while Linker and Dieterich (1992) suddenly changed the normal stress keeping a constant sliding velocity using a servo-controlled system. It is noted that when α = 0, the shear stress reaches a steady-state value instantaneously at the time of the normal stress change. This means that when α is smaller, the shear stress is more sensitive to the normal stress change. On the basis of the assumption that the shear stress is proportional to the true contact area (e.g., Nakatani, 2001), when the normal stress doubles, the contact area becomes twice as large instantaneously if α = 0.

Simulation of earthquake triggering
We use the same model as that by Kato (2004), and Yoshida (2018), who describe the simulation process in detail. A 2D planar fault is loaded at a constant rate V pl in the x-direction, and the fault plane is divided into 256 × 256 square cells with areas of 0.5 km × 0.5 km. The shear stress τ on cell (i, j) is given by: where u(k, l) is the slip at cell (k, l) in the x-direction, K(ik, j-l) is the static shear stress at the center of cell (i, j) due to the uniform unit slip in the x-direction over cell (k, l), V = du/dt is the slip velocity, G is the rigidity, and β is the S wave speed. In all of our simulations, we assumed that V pl = 0.085 m/year, G = 30 GPa, and β = 3.5 km/s. The second term on the right-hand side of Eq. (15) represents the approximate reduction in shear stress due to wave radiation, as introduced by Rice (1993).
We first simulate the unperturbed earthquake cycle using constitutive law (1), normal stress-dependent Nagata law (5), and Eq. (15). Assuming a = 0.03, ab = − 0.012 and 0.0014 (inside and outside the asperity of radius R = 3 km, respectively), σ = 50 MPa, L = 0.003 m, α = 4.2, and c = 20, we obtained the earthquake cycle shown in Fig. 2. Earthquakes occur repeatedly with a recurrence interval T r = 20 years. Bhattacharya et al. (2015) estimated the value of c at 10 to 100 based on laboratory experiments using a rock interface with gouge, which is larger than the value of 2 estimated by Nagata et al. (2012) from bare rock experiments. We obtain α = 4.2 using the scaling relationship with α = 0.2 for the aging law and c = 20. In this section, we apply Coulomb stress perturbations of various  Fig. 2d. Next, we apply dynamic shear and/or normal stress perturbation: The corresponding dynamic Coulomb stress perturbation CFF d (t) is with where S d and N d are the amplitude of the shear and normal stress perturbation, respectively, T p is the period and t = 0 is taken to be the time when the perturbation is applied. Figure 3a shows the Coulomb stress perturbation applied at 7/8 of the recurrence interval T r ~ 20 years, as indicated with the arrow in Fig. 2a. The assumed parameters are S d = 0 MPa, N d = − 1.48 MPa, μ = 0.7, α = 4.2, c = 20, and T p = 100 s. Figure 3b, c shows the slip velocity and displacement, respectively, at the three points in the asperity shown in Fig. 2d. Dieterich (1994) represented the time to instability as a function of slip velocity using a single springblock model and pointed out that an increase in the slip velocity due to a static stress change could trigger an earthquake. This model can be extended to incorporate dynamic stress perturbation, which reduces the frictional strength, causing an increase in slip velocity. Figure 3b shows that the slip velocity after the perturbation ceased is higher than that before the perturbation starts. Figure 4 shows that the slip accelerates after the perturbation starts at t = 0, with the triggered event occurring approximately 17.4 days later. The increase in slip velocity promotes the rupture nucleation growth, resulting in the reduction of the time to instability from 2.5 year to 17.4 days. Figure 5a shows the relationship between the amplitudes of the stress perturbation and the time to instability. The values of the friction parameters are the (16) same as those in Fig. 2. The red circles show the effects of the shear stress change without normal stress perturbation (N d = 0), and the other circles represent the effects of the normal stress changes with S d = 0. The assumed values of α are 0 (orange), 2.1 (green), and 4.2 (black), respectively. The time to instability decreases, i.e., the triggering potential is high, as the stress amplitude increases. As small values of α results in large values of ΔCFF d for negative N d as denoted by Eq. (18), the values of N d required to cause a certain time to instability decrease when α is small. The relationships of ΔCFF d calculated using Eq. (18) and the time to instability are plotted in Fig. 5b, which shows that the curves for the four cases closely overlap. This means that the triggering potential is determined mainly by a value of ΔCFF d . Triggering by static stress change is also simulated. At t = 0, static stress jump is applied: The corresponding Coulomb stress jump is written as Figure 6a shows the effect of the static stress changes S s and N s on the time to instability by assuming μ = 0.7, α = 4.2, c = 20, and the same friction parameters as in Fig. 2a. Figure 6b shows the relationship between ΔCFF s and the time to instability. The two curves for the shear and normal stress jumps completely overlap. In the static case the triggering potential is a function of ΔCFF s . Dieterich et al. (2000) pointed out that a modified Coulomb stress change based on conventional RSF for static triggering is given by ΔCFF m = S s − (μ − α) N s . If scaling relationship (13) is satisfied, ΔCFF m takes the same value as ΔCFF S evaluated using Eq. (21) for given S s and N s .
The change Δ s ln(V/V*) due to the static jump of CFF can be approximately obtained using Eq. (7), . Fig. 4 a Changes in the shear stresses at the three points denoted in Fig. 2d including the period during the triggering perturbation (shown by the expanded time axes in Fig. 3), the following period, and the triggered instability. b Changes in velocities. c Changes in displacements. The parameters are the same as those in Fig. 2. Figure 3 shows expansions of these plots around t = 0 Note that in the previous section, we derived some equations under the condition of V = const. to reproduce the experimental results obtained by Linker and Dieterich (1992). As V is not controlled in the earthquake cycle simulations, lnV increases in proportional to ΔCFF S .
Comparing Figs. 5b and 6b, we find that values of ΔCFF s for static change smaller than ΔCFF d for dynamic change can trigger instability. This means that a value of ΔCFF itself cannot be used for evaluating the triggering potential to compare the static and dynamic perturbations. This is consistent with previous observations. For instance, Kilb et al. (2000) reported that the aftershocks of the 1992 Landers earthquake (M7.3) were possibly triggered not only by static stress changes, but also by dynamic stress changes with much larger amplitude than that of the static changes. In the case of dynamic triggering, the applied Coulomb stress becomes zero after the perturbation ceases, as implied by Eq. (19), and ΔCFF d makes no direct contribution to the first term of the right-hand side of Eq. (7) for changing Δln(V/V*).
Next, we consider the increases of V when the following dynamic perturbation CFF d (t): The velocity when the perturbation is applied, based on Eq. (7), is approximately or where V i is the velocity just before applying the perturbation. The displacement during the perturbation is given by As CFF d (t) = 0 after the perturbation ceases, Eq. (8) approximately leads to Here, we ignored the integral of 1/θ in Eq. (8) because the values were less than 5% of values estimated using the right-hand side of Eq. (26) in the numerical simulations in Fig. 5. Moreover, Δτ has a small value after the perturbation because of the displacement. Heterogeneous displacement inside the asperity can cause both stress release and stress concentration depending on location and heterogeneity.
We evaluate the equivalent ΔCFF eq for the dynamic perturbation, which results in the same velocity increase as that obtained in the case of static stress perturbation when ΔCFF s = ΔCFF eq . Comparing Eqs. (22) and (26), ΔCFF eq is represented by ΔCFF eq is independent of a value of c including c = 0 (aging law) when the scaling relations (13) and (14) are satisfied. Figure 7 shows ΔCFF eq using the average velocity over the asperity before the perturbation starts as V i . It is found that the triggering potentials for the dynamic and static perturbations are similar when ΔCFF s = ΔCFF eq . Effects of stress perturbation timing were briefly discussed by Yoshida (2018) assuming constant normal (23) shear stress Ss normal stress Ns a b Fig. 6 a Relationship of the static stress perturbation amplitude and the time to instability: only shear stress changes are assumed (blue circles), and only the normal stress changes are assumed (the purple squares). b Relationship between ΔCFF s and the time to instability stress. As the effects of timing may not be changed essentially by including the normal stress dependence, we do not discuss this subject here.

CFF due to the 2016 Mie offshore earthquakes
In this section, we will calculate ΔCFF eq using Eq. (27) for the Mie offshore earthquake which occurred around the Nankai Trough on April 1, 2016. Takemura et al. (2018) obtained the moment tensor of this earthquake by inversion analysis using the broadband records from the full-range seismograph network (F-net) stations. The inversion solution was dip = 20º, strike = 250º, rake = 110º, and M w = 5.9. Maeda et al. (2017) developed an open-source integrated parallel simulation code (OpenSWPC) for modeling seismic wave propagation in 3D heterogeneous viscoelastic media. In this paper, we updated OpenSWPC to calculate a stress tensor. We obtain the stress field for the 2016 Mie offshore earthquake using the modified code, assuming a point source with the moment tensor obtained by Takemura et al. (2018) and Japan Integrated Velocity Structure Model (Koketsu et al. 2012).
We calculate Coulomb stresses at grid points with intervals of 10 km on the interface of oceanic layer 2 (basement) of the Philippine Sea Plate. The point adjacent to the epicenter of the Mie offshore earthquake, which is denoted by a star symbol in Fig. 8, is excluded because the point source approximation is not applicable. The total number of the grid points is 15 × 13-1. The point adjacent to the epicenter of the Mie offshore earthquake, which is denoted by a star symbol in Fig. 8, is excluded because the point source approximation is not applicable. The assumed mechanism of the target fault is dip = 15º, strike = 240º, and rake = 110º, which are the same as those Satake (1993) assumed for the 1944 Tonankai earthquake. Figure 9a, b shows shear stress τ d (t) and normal stress σ d (t), respectively, calculated at point A shown in Fig. 8c. Nagata et al. (2012) obtained c = 2 through a laboratory experiment, and Linker and Dieterich (1992) obtained α = 0.2 for the aging law. Using these values, the scaling relationship (13) leads to α = 0.6 for the Nagata law. Figure 9 (c) shows CFF d (t) based on Eq. (23) by assuming these values and μ = 0.6. We calculated CFF d (t) from the occurrence time t = 0 to 450 s (the latter part is not shown in Fig. 9), and we obtained ΔCFF s by averaging CFF d (t) from 400 s to 450 s when the most of oscillations vanishes. Figure 8a, b shows the distributions of the maximum value CFF d_max of CFF d (t) and ΔCFF s , respectively. CFF d_max and ΔCFF s have large values near the hypocenter with largest values of 310 kPa and 38 kPa, respectively, at point A shown in Fig. 8c.
Next, we obtain ΔCFF eq from CFF d (t) using Eq. (27), assuming values of many parameters that are not constrained by observations. Hence, we provide an example of evaluation in this paper. The assumed variation with depth, γ = a − b, is shown in Fig. 10. In the shallow and deep regions, γ 0 = 5.0 × 10 −4 , and in the seismic region, γ seis = − 3.0 × 10 −4 . The variation of γ is similar to that in model 4 assumed by Kato and Hirasawa (1999). They performed a numerical simulation of a large earthquake cycle in the Tokai district using the aging law. The other friction parameters for the aging law are assumed to be a = 5.2 × 10 −4 , L = 1.0 × 10 −2 m, and μ = 0.6, and they are independent of depth. The parameters for the Nagata law are determined using the scaling relation (13) with c = 2 from the aging law parameters.
When R OP is assumed to be 68%, the maximum value of ΔCFF eq at point B is comparable to the maximum value of ΔCFF s at point A. Figure 8c, d shows the distribution of ΔCFF eq and the effective normal stress, respectively. The ΔCFF eq distribution shows a different pattern than the CFF d_max distribution. For a larger R OP , e.g., for a larger pore pressure, ΔCFF eq becomes larger. As Fig. 11  The hypocenters of the low-frequency tremors which occurred from April 1 to April 13 determined by Araki et al. (2017). A SSE associated with these tremors was considered to be triggered by the Mie offshore earthquake. Blue contours show the depth of the plate boundary shows ΔCFF eq as a function of R OP at point B, ΔCFF eq is very sensitive to R OP . This possibly explains why triggered earthquakes have been frequently observed in geothermally active or volcanic areas (e.g., Hill et al. 1993;Aiken and Peng 2014). If the pore pressure is locally very high, earthquake triggering can easily occur in such areas. Wei et al. (2018) also pointed out that low effective normal stress on the shallow subduction interface is needed to reproduce the repeating SSEs.

Discussion and summary
After the 2016 Mie offshore earthquake, shallow SSEs, shallow very low-frequency events, and shallow tremor were observed near the Nankai trough, which were probably triggered by the Mie offshore earthquake (Annoura et al. 2017;Araki et al. 2017;Nakano et al. 2018). Araki et al. (2017) reported SSEs and associated tremor swarms using data from an array of borehole and seafloor instruments. Figure 8e shows the hypocenters of the low-frequency tremors which occurred from April 1 to April 13  Relationship between ΔCFF eq and over-pressure rate R op at point B denoted in Fig. 9c determined by Araki et al. (2017). They considered that a SSE accompanied by these tremors was triggered by the 2016 Mie offshore earthquake. As ΔCFF eq due to the Mie offshore earthquake have large values near the trough, it may have triggered the shallow SSEs. Although the distribution of ΔCFF eq shown in Fig. 8c depends on uncertain parameter values, the large ΔCFF eq value near the trough was mainly caused by low effective normal stress. The shallow region was assumed to have positive a-b in the previous section. Many investigators (e.g., Du et al. 2003;Shibazaki et al. 2012;Wei et al. 2015) performed numerical simulation of SSEs assuming a-b < 0 for the patch at the unstable-stable transient zone where SSEs occur, although Perfettini and Ampuero (2008) indicated that SSEs could occur even if a-b > 0. If we assume a larger b resulting in a-b < 0 in the shallow region, ΔCFF eq would take a larger value in the shallow region as expected from Eq. (27). We confirmed that the ΔCFF eq distribution in the shallow region is similar to that in Fig. 8c by calculating ΔCFF eq assuming a-b = γ seis in the shallow region and R OP = 62%, which is lower than the assumed value in Fig. 8c.
We proposed a normal-stress-dependent Nagata law incorporating the experimental results conducted by Linker and Dieterich (1992). Moreover, we obtained the scaling relationship for α between the aging law and the modified Nagata law was obtained (Eq. (13)). When the perturbation of normal and/or shear stress is applied, dynamic and static ΔCFF for the modified Nagata law can be calculated using Eqs. (18) or (21), respectively.
We performed numerical simulations of earthquake triggering due to static and dynamic perturbations in CFF. The triggering potentials for the static and dynamic perturbation cannot be evaluated by simply comparing CFF values. Even if the value of ΔCFF s is smaller than the maximum value of CFF d (t), the triggering potential of the static perturbation could be higher. To compare the triggering potential of static and dynamic CFF perturbations, we derived ΔCFF eq (Eq. (27)) for the dynamic case. If ΔCFF s = ΔCFF eq , both the static and dynamic perturbations have approximately the same triggering potential. On the basis of the numerical simulations, we confirmed that ΔCFF eq is a good measure of the triggering potential for dynamic perturbation.
We calculated CFF d (t), ΔCFF s , and ΔCFF eq for the 2016 Mie offshore earthquake using revised Open-SWPC. ΔCFF eq along the trench has large values, which may have caused the SSEs which occurred after the Mie earthquake. ΔCFF eq is sensitive to the over-pressure ratio of the pore fluid. It is possible that earthquake triggering easily occurs in areas with locally high pore pressure.