A study of earthquake recurrence based on a one-body spring-slider model in the presence of thermal-pressurized slip-weakening friction and viscosity

Earthquake recurrence is studied from the temporal variation in slip through numerical simulations based on the normalized form of equation of motion of a one-body spring-slider model with thermal-pressurized slip-weakening friction and viscosity. The wear process, whose effect is included in the friction law, is also taken into account in this study. The main parameters are the normalized characteristic displacement, Uc, of the friction law and the normalized damping coefficient (to represent viscosity), η. TR,D, and τD are the recurrence time of events, the final slip of an event, and the duration time of an event, respectively. Simulation results show that TR increases when Uc decreases or η increases, D and τD decrease with increasing η, and τD increases with Uc. The timeand slip-predictable model can describe the temporal variation in cumulative slip. When the wear process is considered, the thickness of slip zone, h, which depends on the cumulated slip, S(t)= ∑ D(t), i.e., h(t)=CS(t) (C is a dimensionless increasing rate of h with S), is an important parameter influencing TR and D. Uc is a function of h and thus depends on cumulated normalized slip, ∑ U , with an increasing rate of C. In the computational time period, the wear process influences the recurrence of events and such an effect increases with C when C > 0.0001. When viscosity is present, the effect due to wear process becomes stronger. Both TR and D decrease when the fault becomes more mature, thus suggesting that it is more difficult to produce large earthquakes along a fault when it becomes more mature. Neither the time-predictable nor the slip-predictable model can describe the temporal variation in cumulative slip of earthquakes under the wear process with large C.


Introduction
Earthquake recurrence that is relevant to the physics of faulting is an important factor in seismic hazard assessment. It is related to the seismic cycle, which represents the occurrence of several earthquakes in the same segment of a fault during a time period. Figure 1 exhibits the general pattern of time variation in slip and particle velocity during a seismic cycle. In the figure, T R is the recurrence (also denoted by repeat or inter-event) time of two events in a seismic cycle, τ D is the duration time of slip of an event, D is the final slip of an event, and V m is the peak value of particle velocity of an event. The four parameters could be constants in the time history when all model parameters do not vary with time and could also vary with time, represented by T R (t), τ D (t), D(t), and V m (t), when one of model parameters does vary with time. Sykes and Quittmeyer (1981) pointed out that the major factors in controlling T R are the plate moving speed and the geometry of the rupture zone. Based on Reid's elastic rebound theory (Reid, 1910), Schwartz and Coppersmith (1984) assumed that an earthquake occurs when the tectonic shear stress on a fault is higher than a critical level, which is dependent on the physical conditions of the fault and the loading by regional tectonics. Since in their work a fault has a homogeneous distribution of physical properties under constant tectonic loading, earthquakes could happen regularly.
Some observations exhibit periodicity for different size earthquakes. Bakun and McEvilly (1979) obtained T R ≈ 23 ± 9 years for M ≈ 6 earthquakes at the Parkfield segment of the San Andreas fault, USA since 1857. Sykes and Menke (2006) Figure 1. A general pattern of time variations in slip (thick solid line) and particle velocity (thin solid curve) during a seismic cycle: T R is the recurrence time or the inter-event time of two events in a seismic cycle; τ D is the duration time of slip of an event; D is the final slip of an event; and V m is the peak particle velocity of an event.
quakes in the Nankaido segments of the Nankai Trough, Japan. Okada et al. (2003) gained T R = 5.5 ± 0.7 years for earthquakes with M = 4.8 ± 0.1 off Kamaishi, Japan, since 1957. Nadeau and Johnson (1998) inferred an empirical relation between T R and seismic moment, M o : To make this relation valid, the stress drop, σ , or the long-term slip velocity of a fault, v l , must be in terms of M o . Based on three data set from eastern Taiwan, Parkfield, USA, and northeastern Japan, Chen et al. (2007) inferred T R ∼ M 0.61 o . Beeler et al. (2001) proposed a theoretical relation: T R = σ 2/3 M 1/3 o /1.81µv l , where µ is the rigidity of the fault-zone materials.
However, the main factors in influencing earthquake occurrences commonly are spatially heterogeneous and also vary with time. Thus, the recurrence times of earthquakes, especially large events, are not constant inferred either from observations (Ando, 1975;Sieh, 1981;Kanamori and Allen, 1986;Wang and Kuo, 1998;Wang, 2005;Sieh et al., 2008) or from modeling (Wang, 1995(Wang, , 1996Ward, 1996Ward, , 2000Wang and Hwang, 2001). Kanamori and Allen (1986) observed that faults with longer T R are stronger than those with shorter T R . Davies et al. (1989) proposed that the longer it has been since the last earthquake, the longer the expected time until the next. Wang and Kuo (1998) observed that for M ≥ 7 earthquakes in Taiwan T R strongly follows the Poissonian processes. Enescu et al. (2008) found that the probability density distribution of T R can be described by an exponential function. From the estimated values of T R of earthquakes happened on the Chelungpu fault in central Taiwan from trenching data, Wang (2005) found that the earthquakes occurred non-periodically.
In order to interpret earthquake recurrences, Shimazaki and Nakata (1980) proposed three simple phenomenological models. Each model has a constantly increasing tectonic stress that is controlled by a critical stress level, σ c , for failure and a base stress level, σ b . The three models are (1) the perfectly periodic model (with constant σ c , σ b , and σ ); (2) the time-predictable model (with constant σ c , variable σ b , and variable σ ); and (3) the slip-predictable model (with variable σ c , constant σ b , and variable σ ). For the first model, both T R and D of next earthquake can be predicted from the values of T R or D of previous ones. For the second model, T R of next earthquake can be predicted from the values of D of previous ones. For the third model, D of next earthquake can be predicted from the values of T R of previous ones. However, debates about the three models have been made for a long time. Some examples are given below. Ando (1975) suggested that the second model worked for post-1707 events but not for pre-1707 ones in the Nankai Trough, Japan. Wang (2005) assumed that the second model could describe the earthquakes occurred on the Chelungpu fault, Taiwan, in the past 1900 years. For the Parkfield earthquake sequence, Bakun and McEvilly (1984) took different models, while Murray and Segall (2002) considered the failure of the second model. From laboratory results, Rubinstein et al. (2012) assumed the failure of the time-and slippredictable models for earthquakes. Some models, for instance the crack model and dynamical spring-slider model, have been developed for fault dynamics, even though the seismologists have not a comprehensive model. There are several factors in controlling fault dynamics and earthquake ruptures (see Bizzarri, 2009;Wang, 2017b). Among the factors, friction (Nur, 1978;Belardinelli and Belardinelli, 1996) and viscosity (Jeffreys, 1942;Spray, 1993;Wang, 2007) are two significant ones.
Some results concerning earthquake recurrence are simply explained below. Erickson et al. (2008) suggested that aperiodicity in earthquake dynamics is due to either the nonlinear friction law (Huang and Turcotte, 1990) or the heterogeneous stress distribution (Lapusta and Rice, 2003). Wang and Hwang (2001) emphasized the importance of heterogeneous frictional strengths. Mitsui and Hirahara (2009) pointed out the effect of thermal pressurization. Dragoni and Piombo (2011) found that variable strain rate causes aperiodicity of earthquakes. Bizzarri and Crupi (2014) found that T R is dependent on the loading rate, effective normal stress, and characteristic distance of the rate-and state-dependent friction law.
As mentioned previously, numerous studies have been made for exploring the frictional effect on earthquake recurrence. However, studies concerning the viscous effect on earthquake recurrence are rare. In the following, we will investigate the effects of slip-weakening friction due to thermal-pressurization and viscosity on earthquake recurrence based on the one-body spring-slider model.
2 One-body model Figure 2 displays the one-body spring-slider model. In the model, m, K, σ n , F , η, u, v (= du/dt), v p , and u o = v p t denote, respectively, the mass of the slider, the stiffness (or spring constant) of the leaf spring, the normal stress, the frictional force between the slider and the moving plate, the damping coefficient (to represent viscosity as explained below), the displacement of the slider, the velocity of the slider, the plate moving speed, and the equilibrium location of the slider. The frictional force F (with the static value of F o ) is usually a function of u or v. Viscosity results in the viscous force, , between the slider and the moving plate, and is in terms of v. A driving force, Kv p t, caused by the moving plate through the leaf spring pulls the slider to move. The equation of motion of the model is When Kv p t ≥ F o , F changes from static frictional force to dynamic one and thus makes the slider move. Among the physical models to approach earthquake faults, the single spring-slider model, which can represent a single fault, is actually the simplest one. However, based on this simple model in the presence of thermal-pressurized friction and viscosity we can obtain good simulations of earthquake recurrences along a single fault. Results can exhibit the frictional and viscous effects on earthquake recurrence. The frictional force F (u, v) is controlled by several factors (see Wang, 2016, and cited references therein). An effect combined from temperature and fluids in a fault zone can result in thermal pressurization (abbreviated as TP hereafter) which would yield a shear stress (resistance) on the fault plane (Sibson, 1973;Lachenbruch, 1980;Rice, 2006;Wang, 2009Wang, , 2011Wang, , 2016Wang, , 2017aBizzarri, 2009). Rice (2006) proposed two end-member models of TP, i.e., the adiabatic undrained deformation (AUD) model and slip-on-aplane (SOP) model. Since the characteristic distance of the Figure 2. One-body spring-slider model. In the figure, t, m, K, η, V p , σ n , F , u, and u o denote, respectively, the time, the mass of the slider, the spring constant, the damping coefficient, the driving velocity, the normal stress, the frictional force, displacement of the slider, and the equilibrium location of the slider (after Wang, 2016). SOP model cannot be associated with the wear process, the SOP model is not used in this study. The AUD model is related to a homogeneous simple strain ε at a constant normal stress σ n on a spatial scale of the sheared layer. Its shear stress-slip function, (Rice, 2006), which decreases exponentially with increas- h, µ f , and are, respectively, the fluid density, heat capacity (in J • C −1 kg −1 ), the thickness, frictional strength, and the undrained pressurization factor of the fault zone. The parameter is , where β f is isothermal compressibility of the pore fluid, β n is isothermal compressibility of the pore space, λ f is isobaric volumetric thermal expansion coefficient for the pore fluid, and λ n is isobaric volumetric thermal expansion coefficient for the pore space.
Based on the AUD model, Wang (2009Wang ( , 2016, 2017a-c) took a simplified slip-weakening friction law (denoted by the TP law hereafter): (2) The value of F (u) at u = 0 is F o , i.e., the static friction force. An example of the plot of F (u) versus u for five values of u c , i.e., 0.1, 0.3, 0.5, 0.7, and 0.9 m when F o = 1N m −2 , which are taken from Wang (2016), is displayed in Fig. 3. F (u) decreases with increasing u and its decreasing rate, γ , decreases with increasing u c . The force drop is lower for larger u c than for smaller u c for the same final slip. When u u c , exp(−u/u c ) ≈ 1 − u/u c , thus indicating that u −1 c is almost γ at small u. This TP law is used in this study.
A detailed description about viscosity and the viscous force (v) can be found in Wang (2016), and only a brief explanation is given below. Jeffreys (1942) first and then numerous authors (Byerlee, 1968;Turcotte and Schubert, 1982;Scholz, 1990;Rice et al., 2001;Wang, 2016) emphasized the viscous effect on faulting due to frictional melts. The viscosity coefficient, υ, of rocks is influenced by T (see Turcotte and Schubert, 1982;Wang, 2011). Spray (2005, and cited references therein) observed a decrease in υ with increas- ing T . He also stressed that frictional melts with low υ could produce a large volume of melting, thus reducing the effective normal stress. This behaves like fault lubricants during seismic slip.
The physical models of viscosity can be found in several articles (e.g., Cohen, 1979;Hudson, 1980). The stress-strain relationship is σ = Eε, where σ and E are, respectively, the stress and the elastic modulus for an elastic body; and σ = υ(dε/dt), where υ is the viscosity coefficient for a viscous body. Two simple models with a viscous damper and an elastic spring are often used to describe the viscous materials. A viscous damper and an elastic spring are connected in series, leading to the Maxwell model, and in parallel, resulting in the Kelvin-Voigt model (or the Voigt model). According to Hudson (1980), Wang (2016) proposed that the latter is more suitable than the former for seismological problems and thus the Kelvin-Voigt model, whose constitution law is σ (t) = Eε(t) + υdε(t)/dt, is taken here and displayed in Fig. 2. The viscous stress is υv.
In order to investigate the viscous effect in a dynamical system, Wang (2016) defined the damping coefficient, η, based on the phenomenon that an oscillating body damps in viscous fluids. According to Stokes' law, η = 6π Rυ for a sphere of radius R in a viscous fluid of υ (see Kittel et al., 1968). Hence, the viscous force in Eq. (1) is represented by = ηv. Note that the unit of η is N(m/s) −1 . Since υ decreases with increasing T , η decreases with increasing T . Hence, η can vary with time during faulting. This point has been studied by Wang (2017b) for the generation of nuclear phase before an earthquake ruptures. In this study, constant η is considered for each case. Some authors (Knopoff et al., 1973;Cohen, 1979;Rice, 1993;Xu and Knopoff, 1994;Knopoff and Ni, 2001;Bizzarri, 2012a;Dragoni and Santini, 2015) considered that viscosity plays a role on causing seismic radiation to release strain energy during faulting.

Normalization of equation of motion
Putting Eq. (2) and = ηv into Eq. (1) leads to Equation (3) is normalized for easy numerical computations based on the normalization parameters, which is dimensionless: The normalized velocity, acceleration, and driving velocity are , respectively. = ω/ω o is defined as the dimensionless angular frequency, and thus the phase ωt becomes τ . For the purpose of simplification, η/(mK) 1/2 is denoted by η below. Substituting all normalization parameters into Eq. (3) leads to In order to numerically solve Eq. (4), we define two new parameters, i.e., y 1 = U and y 2 = dU/dτ . Equation (4) can be rewritten as two first-order differential equations: We can numerically solve Eq. (5) by using the fourth-order Runge-Kutta method (Press et al., 1986). In general, the values of D o are several meters and ω o are in the range of 0.1 Hz to few Hz (see Wang, 2016). This leads to that D o ω o has an order of magnitude of 1 m s −1 . The value of V p must be much smaller than 1 because of v p ≈ 10 −10 m s −1 . Since the value of V p mainly influences the recurrence time, T R , between two events and can only make a very small influence on the pattern of time variations in velocities and displacements of events. In order to study long-term earthquake recurrence, there must be numerous modeled events with clear and visualized time functions of displacements and velocities for an event in the computational time period. If V p = 10 −10 is considered, T R is very long and thus τ D is much shorter than T R . This makes the time function of an event displayed in the long-term temporal variation in slip looks like just a step function for the displacements and an impulse for the velocities. Hence, in order to get fine visualization a larger value of V p is necessary. The value of V p τ is usually very small during an event and cannot influence the rupture because of a very tiny value of V p . Numerical test shows that when V p > 10 −2 , the value of V p τ is not small during an event and can influence the rupture. Hence, V p = 10 −2 is taken in this study.
The backward slip is not allowed in the simulations because of common behavior of forward earthquake ruptures. A phase portrait, which is a plot of a physical quantity, y, versus another, x, i.e., y = f (x), is commonly used to represent nonlinear behavior of a dynamical system (Thompson and Stewart, 1986). The intersection point between f (x) and the bisection line of y = x is defined as the fixed point, that is, f (x) = x. If f (x) is continuously differentiable in an open domain near a fixed point x f and |f (x f )| < 1, attraction can appear at the fixed point. Chaos can also be generated at some attractors. The details can be seen in Thompson and Stewart (1986). In this study, the phase portrait is the plot of V /V max versus U/U max .

Simulation results
Numerical simulations lead to the temporal variations in particle velocities and displacements as displayed in Fig. 1 Simulation results displayed in these figures show that the maximum values of both V and U decrease from case (a) to (d) in each figure. Hence, the maximum velocity and maximum displacement, which are denoted by V max and U max , respectively, for case (a) can be taken as the scaled factor to normalize the waveforms from case (a) to (d). This makes it easy to compare the waveforms of the four cases in each figure.
The cases excluding the viscous effect, i.e., η = 0, are first simulated and results are shown in Fig. 4 for four values of  Fig. 6, and 0.8 in Fig. 7.
The left panels in Fig. 4 with η = 0 show that the peak velocity of an event, V m , and final slip, D, with the re-spective maximum values in case (a) as mentioned above, for all simulated events decrease with increasing U c . From Fig. 3, the force drop, F , decreases with increasing U c for a certain final slip, thus indicating that larger F yields higher V m and larger D. This interprets the negative dependence of V m and D on U c . When the viscous effect is absent, i.e., η = 0, the value of τ D increases with U c , while T R decreases with increasing U c . When U c = 1, V m and D are both very small and the system behaves like creeping of a fault. In the right panels, there are two fixed points for each case: one is called the non-zero fixed point at larger V and larger U and the other the zero fixed point at V = 0 and U = 0. The slope at a fixed point is defined to be d(V /V max )/d(U/U max ) = (U max /V max )(dV /dU ). The absolute values of slope at the two fixed points decrease with increasing U c , thus suggesting that the fixed point is not an attractor for small U c and could be an attractor for larger U c . The phase portrait for U c = 1 is very tiny, because the final slip for U c = 1.0 is much smaller than those for U c = 0.2, 0.4, and 0.8. Hence, U c = 1 will not be taken into account in the following simulations.
The left panels in Figs. 5-7 show that V m and D decrease when either U c or η increases, while τ D increases with η and U c . Meanwhile, T R increases when either η increases or U c decreases. The right panel shows that the phase portraits coincide for all simulated events for a certain η. The absolute values of slope at the two fixed points decrease when either U c or η increases. This suggests that the fixed point is not an attractor for small U c and low η, and can be an attractor for large U c and high η. Like Fig. 4, the final slip decreases with increasing U c .
In Figs. 5-7 we can see that the temporal variation in cumulative slip can be described by the perfectly periodic model as mentioned above. Hence, when U c and η do not change with time, the rate of cumulative slip retains a constant in the computational time period. This is similar to the simulation results with the periodical earthquake occurrences obtained by some authors (e.g., Rice and Tse, 1986;Ryabov and Ito, 2001;Erickson et al., 2008;Mitsui and Hirahara, 2009) based on the one-body model with rate-and statedependent friction or velocity-weakening friction. However, the present result is inconsistent with the simulation results, from which either the time-predictable model or the slippredictable model cannot interpret the temporal variation in cumulative slip, based on the same model obtained by others (e.g., He et al., 2003;Bazzarri, 2012b;Bizzarri and Crupi, 2014;Kostić et al., 2013a, b;Franović et al., 2016). The differences between the two groups of researchers might be due to distinct additional constrains in respective studies. Although the detailed discussion of such differences is important and significant, it is out of the scope of this study and ignored here.
The phase portraits in Figs. 5-7 exhibit two kinds of fixed points as mentioned above. The absolute values of slope at the non-zero fixed point are higher than 1 and decreases with  increasing η. This means that larger η is easier to generate an attractor than small η. However, the reducing rate of absolute value of slope decreases with increasing U c . The absolute values of slope of the zero fixed point are higher than 1 and decrease with increasing η. This suggests that the zero fixed points can be an attractor. This behavior becomes weaker when U c increases. Figures 4-7 show that when U c and η are constants during the computational time periods, the general patterns of temporal variations in cumulated slip do not change. Some of the previous studies (e.g., Bizzarri, 2012a, b;Franović et al., 2016) suggest that the patterns of temporal variations in cumulated slip can change with time. The changes of U c and η with time should play the main roles. From u c = ρ f C V h/µ f of the TP model (see Rice, 2006), the width of the slipping zone, h, where the maximum deformation is concentrated (Bizzarri, 2009), is a significant parameter in this study. From geological surveys, Rathbun and Marone (2010) observed that h is not spatially uniform even within a single fault. Hull (1988) and Marrett and Allmendinger (1990) found that the wear processes occurring during faulting could widen h, and thus h could vary with time. According to the results gained by several authors (e.g., Power et al., 1988;Robertson, 1983;Bizzarri, 2010), Bizzarri (2012b) proposed a linear dependence of h on the cumulated slip, S(t) = D(t), i.e., h(t) = CS(t), where C is a dimensionless increasing rate of h with S and is consid-ered to be a constant in each case. Based on h(t) = CS(t), the more mature the fault is, the thicker its slip zone is. Since u c is proportional to h and U c = u c /D o , U c is related to C. Here, we assume that U c varies with cumulative slip in the following way: U c = U co + C U (t), where U co is the initial value of U c . Simulation results for four values of C are shown in Figs. 8-12: (a) for C = 0.0001, (b) for C = 0.001, (c) for C = 0.01, and (d) for C = 0.05 when η = 0 in Figs. 8-10; (a) for C = 0.0001, (b) for C = 0.001, (c) for C = 0.01, and (d) for C = 0.038 when η = 1 in Fig. 11; and (a) for C = 0.0001, (b) for C = 0.001, (c) for C = 0.01, and (d) for C = 0.0136 when η = 1 in Fig. 12. The initial values of U c are 0.1 for Fig. 8, 0.5 for Fig. 9, 0.9 for Fig. 10, 0.1 for Fig. 11, and 0.5 for Fig. 12. Note that the value U c varies with time due to time-varying h.
The left panels of Figs. 8-12 show that V m , D, τ D , and T R are all similar when C ≤ 0.001. However, in general V m and D decrease with increasing C, T R slightly decreases with increasing C, and τ D slightly increases with C. A decrease in D is particularly remarkable when C ≥ 0.01. When h is wider than a critical value with C = 0.05 for η = 0, normal earthquakes cannot occur and only creeping may happen. The critical value of h decreases when the viscous effect is present with η = 1 in this study. This decrease is also influenced by U c : C = 0.038 when the initial value of U c is 0.1 and C = 0.0136 when the initial value of U c is 0.5. Obviously, T R decreases with increasing C, thus leading to a de-crease in T R with increasing h. This is similar to the result obtained by Bizzarri (2010Bizzarri ( , 2012b. However, the viscous effect was not included in his studies. The right panels of Figs. 8-12 show that the phase portraits for C = 0.001 are slightly different from those for C = 0.0001 even though the patterns of their variations in V and U are similar, while the phase portraits for C > 0.001 are different from those for C ≤ 0.001. An increase in h due to an increase in C with cumulative slip enlarges U c . This can be explained from Fig. 3, which shows that larger U c yields a lower F than smaller U c for the same final slip. Hence, an increase in U c produces a decrease in F , thus resulting in low V m and small D. In addition, an increase in U c makes exp(−U/U c ) approach unity, especially for small U , thus reducing the nonlinear effect caused by TP friction.
Unlike Figs. 4-7, the size of phase portraits in the right panels of Figs. 8-12 decreases with increasing C. This reflects a decrease in both T R and D of events with increasing C as mentioned previously. The absolute values of slope at the non-zero fixed point are higher than 1 and only slightly decrease with time when C < 0.01, while the values remarkably decrease with time when C ≥ 0.01. The absolute values of slope at the zero fixed point are higher than 1 and only slightly decrease with time when C < 0.01, while those decrease with time when C ≥ 0.01. Results suggest that the non-zero fixed points for all cases in the study are not an attractor, and those at the zero fixed points can evolve to an attractor with time when C ≥ 0.01. The phenomenon is particularly remarkable for C = 0.05 in Figs. 8-10, C = 0.0380 in Fig. 11, and C = 0.0136 in Fig. 12, and the evolution is faster for large U c than for small U c .

Discussion
The simulation results as mentioned previously demonstrate that when U c and η are constants during the computational time periods, the general patterns of temporal variations in cumulated slip cannot change. In order to investigate the effect of time-varying η and U c on the patterns of temporal variations in cumulated slip, we must consider changes of U c and η with time. The viscosity coefficient can actually vary immediately before and after the occurrence of an earthquake (see Spray, 1883Spray, , 2005Wang, 2017b, c). However, a lack of long-term variation in η does not allow us to explore its possible effect on the change of general patterns of temporal variations in cumulated slip. Here, only the possible effect due to time-varying U c .
As mentioned above, the equality Obviously, U c is controlled by six factors, i.e., ρ f , C V , h, µ f , D o , and . Since the tectonics of a region is generally stable over a long time, the value of D o = F o /K could not change too much and thus would not influence U c . The Debye law (cf. Reif, 1965) gives C v ∼ (T + 273.16) 3 , where 273.16 is the value to convert temperature from Celsius to kelvin, at low T and C v ≈ constant at high T . The threshold temperature, from which C v begins to approach a constant, is 200-300 • K. In this study, C v is almost a constant because of T > 250 • C = 523.16 • K, which is the average ambient temperature of fault zone with depths ranging from 0 to 20 km. Hence, C v is almost constant during a long time and thus cannot influence U c .
The frictional strength, µ f , is influenced by several factors including humidity, temperature, sliding velocity, strain rate, normal stress, and thermally activated rheology (Marone, 1998;Rice, 2006), and thus can change with time (Sibson, 1992;Rice, 2006). Hirose and Bystricky (2007) observed that serpentine dehydration and subsequent fluid pressurization due to coseismic frictional heating may reduce µ f and thus promote further weakening in a fault zone. The pore fluid pressure exists in wet rocks, yet not in dry rocks. Clearly, the time variation in µ f can affects the earthquake recurrences. However, a lack of long-term observations of µ f during a seismic cycle makes the study of its effect on earthquake recurrence be impossible.
The fluid density ρ f and the porosity n depend on T and p. Although there are numerous studies on such dependence (Lachenbruch, 1980;Bizzarri, 2012b; and cited references therein), observed data and theoretical analyses about the values of ρ f and n during a seismic cycle are rare.
The porosity is associated with the permeability, κ. Bizzarri (2012c) pointed out that the time-varying permeability, κ(t), and porosity of a fault zone (cf. Mitsui and Cocco, 2010;Bizzarri, 2012b) can reduce T R . One of the Kozeny-Carman relations (Costa, 2006, and references cited therein) is κ(t) = κ C ϕ 2 (t)d 3 (t)/[1 − ϕ(t)] 2 , where κ C is a dimensionless constant depending on the material in consideration; ϕ is V voids /V tot , where V voids and V tot are, respectively, the pore volume and the total volume of the porous materials; and d is the (average) diameter of the grains, ranging between 4 × 10 −5 and 1 × 10 −4 m (Niemeijer et al., 2010). Usually, κ, ϕ, and d can vary in the fault zone (Segall and Rice, 1995). After faulting κ and ϕ would change and d becomes smaller because of refining of the grains. According to this relation, Bizzarri (2012b) found that κ(t) could significantly reduce T R in comparison with the base model with constant κ. The reason is explained below. An increase in permeability can result in an increase in pore pressure, p f . This can reduce the frictional resistance from τ = µ(σ n − p f ) and thus could trigger earthquakes earlier. Hence, the time-varying permeability can change T R . Nevertheless, we cannot investigate its influence on earthquake recurrence here because there is a lack of a long-term observation of hydraulic parameters during a seismic cycle.
It is significant to explore the factors that can yield a nonperfectly periodic seismic cycle. The width of the slipping zone, h, can be a candidate as pointed out by some authors (e.g., Bizzarri, 2009;Rathbun and Marone, 2010). Since the displacement along a fault is controlled by the fault rheol-ogy, h should depend on the rheology on the sliding interface. The wear processes occurring during faulting could widen h (Hull, 1988;Marrett and Allmendinger, 1990). According to the results gained by several authors (e.g., Power et al., 1988;Robertson, 1983;Bizzarri, 2010), simulation results for various values of C and the results are shown in Figs. 8-10 with η = 0 and in Figs. 11 and 12 with η = 1. Results show that when C > 0.0001, the wear process affects the recurrence of slip and the effect increases with C and when C is larger than an upper-bound value, larger-sized events cannot occur and the earthquake recurrence does not exist. Both T R and D decrease when the fault becomes more mature due to a thicker slip zone. Meanwhile, the viscous effect can also play a secondary role on the earthquake recurrence because it makes upper-bound value become smaller. Although either the time-or slip-predictable model can describe the temporal variations of cumulative slip of events occurring in the earlier time period, they cannot interpret those of events in the later parts. This might suggest that it is more difficult to produce large earthquakes along a fault when it becomes more mature, especially for the cases with viscosity. This implicates that seismic hazard is higher for a young fault than a mature one. Hence, it is significant and important to identify the width of slip zone of an earthquake fault for seismic hazard estimates.

Conclusions
To study the frictional and viscous effects on earthquake recurrence, numerical simulations of the temporal variations in cumulative slip have been conducted based on the normalized equation of a one-body model in the presence of thermal-pressurized slip-weakening friction and viscosity. The wear process, which is included in the friction law, is also taken into account. The model parameters of friction and viscosity are represented, respectively, by U c and η, where U c = u c /D o is the normalized characteristic distance and η is the normalized damping coefficient. Numerical simulation of the time variations in V /V max and cumulative slip U/U max as well as the phase portrait of V /V max versus U/U max are made for various values of U c and η.
Results show that both friction and viscosity remarkably affect earthquake recurrences. The recurrence time, T R , increase when η increases or U c decreases. The final slip, D, and the duration time of slip, τ D , of an event slightly decrease when η or τ D increases and slightly increases with U c . Considering the effect due to wear process, the thickness of slip zone, h, that depends on the cumulated slip, S(t) = D(t), i.e., h(t) = CS(t) (C is a dimensionless constant), is an important factor in influencing earthquake recurrences. U c increases with U with an increasing rate of C. When C > 0.0001, the wear process influences the recurrence of slip and the effect increases with C. When C is larger than an upper-bound value, larger-sized events cannot occur and the earthquake recurrence does not exist. If the slip zone becomes thicker, the fault is more mature. This makes T R and D become shorter. This might suggest that it is more difficult to produce large earthquakes along a fault when it becomes more mature. This phenomenon becomes remarkable when the viscous effect exists because the upper-bound value becomes smaller. The temporal variation in slip cannot be interpreted by the time-predictable or slip-predictable model when the fault is affected by wear process with large C. In addition, the size of phase portrait of V /V max versus U/U max decreases with increasing C. This again reflects decreases in both T R and D of events with increasing C as inferred from the temporal variations in cumulative slip.
Data availability. No data sets were used in this article.
Competing interests. The author declares that he has no conflict of interest.