The variations of characteristic scale length of friction evolution affect earthquake dynamics

Within a fault governing model the characteristic scale length is one of the most relevant physical parameters because it accounts for the so – called fracture energy (density) of the system, its dynamics, the time during which the accumulated stress is released and the seismic waves are excited, the amount of slip developed during an instability event. Friction laboratory experiments reveal that it is not a material property, but that it changes with the sliding velocity. We propose two rather different analytical models to fit laboratory evidence and we incorporate them into a fault model able to simulate repeated earthquakes in the framework of various formulations of rate and state friction. We demonstrate that temporal variations of the scale length do not prevent the system to reach its limit cycle, but they systematically reduce the magnitude of the expected event (both in term of developed slip, and thus seismic moment, and released stress) and also reduce the inter–event time (recurrence interval). Depending on the friction model, the system can penetrate into the stable regime and can either continue the accelerating phase toward to failure or decelerate and abort instability.


Introduction
There is no doubt that modern seismology still seeks improvements in the formulation of models to properly describe the physical processes occurring during faulting, in order to be able to explain the causes of the emission predictions and proposed models against data recorded during seismic sequence and major failure episodes.
Of course, one of the five major pillars illuminating the modern day researches [e.g., Bizzarri, 2014] is represented by laboratory inferences, despite the severe limitations mentioned above. Although crustal earthquakes are thought to be a mixture of fracture and friction phenomena [e.g., Ohnaka, 2003], the prominent part of these laboratory works rely on friction experiments, i.e., on the measure of the rheological properties of preexisting sliding surfaces, in mutual contact and subject to different loading conditions [see sections 2.3 and 2.4 of Bizzarri, 2014 and references therein for a review of this subject]. The pioneering experiments performed by Dieterich [1978] and Ruina [1983] has been followed by several others, culminated in the so-called high velocity rotary shear experiments [see, among others, Sone and Shimamoto, 2009].
All these results highlighted the fact that the pivotal quantity which controls the physics of the system (and therefore the evolution of a natural fault) is represented by the scale length. This quantity -sometime invoked, improperly, as critical distance -basically represents the amount of cumulative slip that the system develops during the breakdown phase, i.e., during the portion of the coseismic stage during which the fault traction degrades from the yield level to the kinetic one. This scale length is reported as characteristic slip-weakening distance (denoted by the symbol d 0 ) in the framework of the slip-dependent governing model, such as the canonical, linear slipweakening model postulated by Andrews [1976], or as the scale length L (or D c as preferred by experimentalists) within the context of the rate-and state-dependent friction laws [Ruina, 1983 and references therein). As clearly stated by Bizzarri and Cocco [2003], d 0 and L are intimately related one to the other, because the description of the weakening process predicted by the two different classes of governing models (slip-dependent or the rate-and state-dependent) are not antithetic, but they can be easily reconciled on theoretical basis [see also Bizzarri, 2011a].
Independently of the assumed friction model, at a general level the scale length represents the possibility to move from the world of the laboratory (where the typical scales of the samples are of the order of tens of cm at maximum) to the world of the natural fault systems (geophysical structures which are expected to be even several tens of km long). Following the findings of Marone and Kilgore [1993] we can assume that the rheological processes acting at a natural scale operate with a larger characteristic length compared to that controlling the laboratory processes. In this framework, a recent study by Leeman et al. [2018], extending previous results by Mair and Marone [1999], states that the characteristic length is not a constant property of the sliding material, but it depends on the sliding speed of the system.
In the present study we propose two different analytical models (see section 2) which fit the laboratory data of Mair and Marone [1999]. By adopting a mass-spring model to solve the elastodynamic problem (see section 3) we incorporate these models in various formulations of rate-and state-dependent friction laws, in order to explore the role of time-variable L in the entire cycle of a natural fault (see sections 4 and 5).

Time variable scale length: the proposed models
As stated above, experimental data by Mair and Marone [1999] clearly indicate that the scale length L is not a constant depending only on the material, acting normal stress and other experimental conditions, but it depends on the sliding speed. Indeed, Figure 9 of Leeman et al. [2018] shows that L systematically increases with velocity and it can vary more than a factor 10 when sliding speed changes from 0.01 mm/s to 10 mm/s. We report the results of Mair and Marone [1999] as full black dots in Figure 1; the inset plot refers to the same velocities considered in their experiments.
In order to fit their data we assume two rather different analytical functions. The first is a logarithmical form [as suggested by Leeman et al., 2018], written as it follows: (1) where L ref is the reference value of the characteristic length (i.e., the value assumed in a model without variations of L), v is the slip velocity, v * is a reference value of the slip velocity and and are the two dimensionless free parameters. Data of Mair and Marone [1999] are fitted very well (with R 2 = 0.95) with = 2 and = *[ln( *) + 11.6].
(Note that and are of the same order of magnitude and therefore they are equally relevant in equation (1)).

Andrea Bizzarri
The other expression we consider is a power law form which reads: (2) where again and are dimensionless free parameters. With = 0.4 and ″105.4 * we also obtain a good fit of data (with R 2 = 0.98); the model (2) is superimposed as blue curve in Figure 1.
Practically, to avoid to have estimates of L lower than L ref (or even negative within model (1)) we consider the two following formulations: ( 3) and (4) respectively. In other words, the assumed models (3) and (4) predict that L = L ref for low slip velocities (namely for v lower than roughly 10 μm/s), while for greater values of the scale length changes with (and therefore implicitly with time on which explicitly depends).
In Figure 2 we superimpose, for models (3) and (4), the expected value of L for different cases of L ref , ranging from 1 mm to 1 cm; these are values commonly expected for real faults [Marone and Kilgore, 1993]. It is clear that, for all cases, for high values of v the power law model (equation (4)) predicts a more relevant change in L compared to the logarithmical model (equation (3)).
represents the speed of a tectonic plate which loads the seismogenic region. The model is represented schematically in Figure 3.
In the framework of this analog 1-D fault model the elastodynamic problem simply is the equation of motion of a damped harmonic oscillator, which can be written as [Bizzarri, 2012c]: where dot indicates the time derivative and (6) U ≡ and F ≡ .  (4)). The inset plots reports L in μm and in a log scale. Note that for v lower than roughly 10 μm/s L is constant (and equal to L ref ) for both the models. For completeness, we also report the values of v c (discriminating between quasi-static and fully dynamic behavior) and v l defining the beginning of the seismic range (see section 3.2 and Table 1 for details).
Mathematically, U in equation (6) indicates the state of the dynamical system, at a generic time t, in the phase space (u, v), u denotes the amount of fault slip at time t and v is its time derivative. Contrarily to the extended fault models, in the present model the discontinuity interface is formally absent, however, within the coseismic phase (defined by the "slip" event), the quantity u can be interpreted as the fault slip, i.e., the discontinuity of the displacement vector in continuum models of faulting [see equation (3) in Bizzarri, 2011a]. On the other hand, the velocity v varies from the very low values of the interseismic phase (or the "stick" phase of the dynamical system, where v ~ v load or less) to the high values attained during the coseismic stage (where v is in so-called seismic range, where v ~ several m/s).
Finally, the traction τ in (6) mathematically represents the fault governing model, which expresses analytically the dependence of the shear components of the stress tensor on the physical observables [u, v, some state variables, etc.; see equation (10) of Bizzarri, 2011a]. In the present paper we will adopt two different friction models, as described in sections 4 and 5.

Numerical approach
As discussed in Bizzarri [2012a], the equation of motion (5) possesses closed-form analytical solutions in some special cases, but it does not in general conditions. Since we will employ non linear governing models (see next sections 4 and 5), it has to be solved numerically. To this goal we use a Runge-Kutta algorithm at fourth order with auto-adaptive time stepping [Press et al., 1992]. We use the algorithms RKQC and RK4 of Press et al. [1992], which has been proved to agree [e.g., Kato, 2001;de Lorenzo and Loddo, 2010] with the results obtained with the ODE45 routine in Matlab ® .

5
Characteristic length and earthquake dynamics The behavior of the dynamical system is characterized by two different regimes, discriminated by the threshold velocity v c = 0.01 mm/s, imposed as a model' s parameter in all the numerical experiments performed here; one is a quasi-static range, for which we neglect inertial effects (i.e., the acceleration is assumed to be negligible, because the fault is evolving very slowly), and the other is the dynamic one, for which the complete problem (equations (5) and (6)) is solved.
We define the occurrence of a dynamic instability (formally the "slip" event or the earthquake) as the time level t [n] where the v exceeds the threshold value of v l = 0.1 m/s, in agreement with previous literature [e.g., Day et al., 2005;Rubin and Ampuero, 2005;Bizzarri and Belardinelli, 2008;Bizzarri and Spudich, 2008;Bizzarri, 2012c]: n being the integer counting the number of instability during the entire seismic cycle. The time interval which separates two subsequent instabilities naturally defines the inter-event time, T cycle (referred also as the seismic cycle): (8)

Stability of the system
The theory predicts a bifurcation in the evolution of the system depending on the dimensionless parameter [e.g., Gu et al., 1984;Rice, 1983], also known as stiffness ratio: when < 1 (i.e., k < k c ) the system is inherently unstable, in that it stores elastic energy during the interseismic loading and then it releases it in a dynamic fashion. On the contrary, when > 1 (i.e., k > k c ) the system slides stably under quasti-static loading, in that the elastic unloading balances the frictional weakening. It can become unstable if the system is subject to a sufficiently strong dynamic loading [e.g., Scholz, 1998]. In the framework of the rate-and state-dependent friction laws the critical stiffness k c , neglecting low order terms, is expressed as: where a and b are two dimensionless governing parameters, which depend on the material and eventually also on temperature, pressure, etc.; the first accounts for the direct effect of velocity on the friction coefficient and the second one controls the releasing phase of stress. More precisely, the difference (B -A) ≡ (b -a) n eff quantifies the decrease of the fault traction from its initial value to the residual level caused by a sudden increase in the loading velocity of a factor e.
Theoretically, the temporal variation of L, as predicted by models (3) and (4), has relevant implication in the degree of instability of the system. Indeed, it can happen that the system, starting from an inherently unstable configuration, enters into the stable behavior, in that the variation of sliding speed causes be greater than unity  (3) and (4), respectively. This means that, when L is not constant, a velocity-weakening regime [defined by the condition b > a; Gu et al., 1984;Rice, 1983], starting with inherently unstable conditions ( < 1), is not a sufficient condition to obtain a sequence of repeated seismic instabilities; we will discuss this result in the next sections. It is interesting to note that while v c and v l are model parameters, assigned a priori before the numerical experiments, v CS depends on the adopted constitutive parameters and on the model adopted to describe the variation of L.

Results with Dieterich-Ruina friction law
Within the framework of Dieterich-Ruina (DR thereinafter) friction the state variable identifies the average contact time of asperities at a microscopic level [Dieterich, 1978;Beeler et al., 1994]: The parameters in the present paper are tabulated in Table 1. Since the adopted model is 1-D, the assumed value of effective normal stress has to be intended as the condition at the earthquake hypocenter. Variations of n eff as such proposed by Linker and Dieterich [1992] and by Dieterich and Linker [1992] are disregarded here, where we assume that the effective normal stress is constant, although they have been considered in more elaborated models [Bizzarri and Cocco, 2006a;2006b;Rice, 2006]. The initial state of the system is ( 0 , 0 ) = ( * , ss (v * )), where the symbol ss stands for the steady state value of friction (defined by the condition (d/dt) = 0, so that ss ( ) = ss ( ) n eff = * -( -) ln ). With the assumed values of v * , and we have: = -20.54 and = 1.67 × 10 -2 .
(a) The system starts at t = 0 from its steady state (at a generic time t * the steady state is defined by the condition ₌ * = 0. slip behavior of the system, which produces repeated instabilities of the same magnitude (Figure 4a), releasing the same amount of stress (Figure 4b), is clearly visible from that figure. The system reaches its limit cycle, because maximum and minimum values of the traction during subsequent instability events are identical; in the correspondence of each single instability the drop in the stress is the same and therefore we can state that, on this fault and in these special conditions, all the occurring earthquakes exhibit the same magnitude. We can also see that the interseismic phase of stress recovery (namely the restrengthening stage) is always the same. This is simply because the state of the system at the end of the deceleration phase is identical over the whole sesimic cycle (see Figure 4c and the inset in Figure 5a); the minimum values of the slip velocity are the same and consequently also the restrengthening phase is the same. In other words, the time required to the fault to completely recover its stress and finally reach the subsequent event is constant. All these results are synthetically represented in the phase portrait ( Figure 4c); at the beginning there is a transient phase, consisting into few instabilities, and then the system enters in its limit cycle, in that the trajectory followed by the system is unique and covered forever. This behavior, which somehow identifies the idealized concept of characteristic earthquake, is expected in these idealized conditions, although it is not a general rule, as discussed by Bizzarri [2012c].
The presence of temporal variations in the characteristic length does not destroy the cyclicity of the instabilities; also with variable L the system undergoes to its limit cycle, but it has severe implications in the whole dynamics of the events. In particular, we can observe for both the models (3) and (4): I. the events are more frequent (in particular for model (4)); the inter-event time T cycle (see equation (8)) is reduced (Figures 4a and 4b). From 170 yr in the reference case, the seismic cycle becomes 147 yr and 94 yr for models (3) and (4), respectively (the variation therefore is of about 13.5 % and 44.7 %, respectively); II. the events have lower magnitude (Figures 4a and 4b), in that the amount of slip developed and the stress released during each instability are lower compared to the reference case. The slip at each instability is 1.7 m for the reference case and it is reduced to 1.4 m and 0.9 m for models (3) and (4), respectively. The stress release (namely, the difference between the peak traction and the residual traction after the total release) is 16.73 MPa in the reference situation and 14.49 MPa and 9.09 MPa in case of models (3) and (4), respectively.
By looking at Figure 5 we can clearly recognize that the evolution of the system is the same at low velocities (as expected, since for both models (3) and (4) L = L ref for v < 10 μm/s, as discussed in section 2). Therefore the time evolutions of v and are identical in this range. Once L starts to vary, the state variable evolution begins to change compared to the reference case (see Figure 5b); this in turn causes a reduction in the peak slip velocity. This is not surprising, since an increase in L reflects into a decrease of the degree of instability of the system (as also predicted by the increase in the stiffness ratio ; see section 3.3). At the same time, the stress release is lower and thus the interseismic phase begins earlier (see Figure 4c). The more relevant change obtained with model (4) compared to that obtained with model (3) can be clearly explained by the fact that in the former case the system spontaneously enters in the stable regime ( > 1), while in the latter case remains in the unstable range (see Figure 5a).
The same conclusions also hold for other values of a and b (nor reported here for brevity) and of L ref . In Figure 6 we For greater values of L ref the system does not enter in its limit cycle, i.e., it is not able to develop instabilities; once the velocity starts to increase (due to the applied external loading) also L starts to increase and this results in an increase of the stiffness ratio , which finally impedes the dynamic acceleration to take place. In other words, with greater values of L ref models (3) and (4) inhibit the occurrence of instabilities.

Results with Ruina-Dieterich model
In this section we adopt a different temporal evolution of the frictional resistance, which follows the Ruina-Dieterich (RD henceforth) law [Ruina, 1983;Beeler et al., 1994;Roy and Marone, 1996]: In this formulation the analytical expression of the friction coefficient is mathematically the same, has again the dimension of time, but it evolves in a different fashion, usually known as slip law (compared to the so-called ageing law of the DR model). This law is also known to predict a behavior of the simulated fault which is globally more unstable than the DR model, for the same adopted constitutive and model parameters and initial conditions. In particular, it predicts a more abrupt stress drop [see also Bizzarri and Cocco, 2003] and a faster dynamics.

This is confirmed by comparing Figures 4 and 7; in the reference cases (black curves) the unperturbed T cycle
is now 138 yr (instead of 170 yr for the DR constitutive model). Correspondingly, the value of traction reached after the instability remains higher in this case and therefore the interseismic restrengthening phase is shorter (compare Figures 4c and 7c).
In contrast to the DR case (Figures 4 and 5), when the spring-slider is governed by the RD governing model it is less affected by the temporal variations of the length scale. Indeed, this occurs because the slider is more unstable than that governed by the DR friction law and because the state variable evolves faster. By introducing a variable L we still observe a reduction in the inter-event time (we now have T cycle = 136 yr for model (3) and T cycle = 134 yr for model (4)), but the variations with respect to the reference case are now only about 1.45 % and 2.9 %, respectively. This is due to the fact that the minimum value of slip velocity (see Figure 7a) is practically unchanged in the three cases and that the traction exhibits a very small variation in its minimum value (see Figure 7b). As a consequence, the interseismic fault restrengthening stage is very similar in all cases and it takes roughly the same times (with the feeble differences of 2 and 4 yr mentioned above).
Things tend to change for increasing values of L ref . When L ref = 1 cm ( Figure 8) the slider, which now starts with ref = 0.37, leaves the quasi-static evolution, enters into the dynamic range (we recall that this occurs for v > v c = 0.01 mm/s; see section 3.2), and -due to the intrinsic high degree of instability guaranteed by the RD model -it enters in the stable range ( > 1) for both models (3) and (4). Indeed, it is able to significantly penetrate the stable range (for v > v CS ; see equations (11) and (12)) and finally it generates seismic instabilities (for v > v l = 0.1 m/s; see section 3.2). We emphasize here that with the DR governing model, for L ref = 2 mm, model (4) does not allow the system to exceed v l (see blue curve in Figure 6a).
With the RD law and L ref = 1 cm the variations of L are so relevant (see also Figure 2) that the slider is significantly affected by them. Now the system continues to reach its limiting cycle, but the inter-event time changes from 125 yr for the reference case to 96 yr and 55 yr for models (3) and (4), respectively (so that the variation is about 23.2 % and 56 %, respectively). Indeed, in this configuration the maximum slip velocity is reduced (especially for model (4); see Figure 8a) and also its evolution in the low velocity range is markedly changed. In the same line, there is a severe reduction of the stress release (see Figure 8b), so that the fault takes less time to recover its strength (see Figure 8c).  (c) Phase portrait. Note that now both models (3) and (4) predict that the system can enter into the stable regime; indeed, in this case we have: v CS = 7.83 x 10 -3 m/s and v CS = 5.87 x 10 -3 m/s (equations (11) and (12), respectively).

Andrea Bizzarri
14 Figure 8. The same as in Figure 7, but now for L ref = 1 cm. In this case both models (3) and (4) cause the system to enter in the stable range, since v CS = 3.54 x 10 -3 m/s and v CS = 1.05 x 10 -3 m/s (equations (11) and (12), respectively).

Discussion and conclusions
We have proposed two analytical models to fit laboratory data performed by Mair and Marone [1999]; see also Leeman et al. [2018]. These friction experiments indicate that the scale length L -which controls the state variable evolution within the context of the rate-and state-dependent governing models [Ruina, 1983] -is not a property which depends only on the type of material, but it also depends on the sliding speed (formally on the fault slip velocity). At a more fundamental level, in a dynamic model of faulting processes, the scale length definitively is one of the most relevant physical parameters, in that it controls the whole dynamics of the system, in terms of time scales of the stress release, amount of slip developed during an instability event, fracture energy density spent to break a new portion of the fault (and thus it also controls the seismic wave radiation).
We have shown (see Figure 1) that experimental data can be fitted both with a logarithmic form (model (1)) and with a power law form (model (2)). We extend these models outside the range of sliding speeds explored in laboratory experiments by postulating models (3) and (4); see Figure 2. We incorporate these proposed models into two different analytical expressions of the rate-and state-dependent friction laws, the Dieterich-Ruina (DR) model (equation (13)) and the Ruina-Dieterich (RD) model (equation (14)). Then we solve the elastodynamic problem for a single, isolated, 1-D spring-slider model (namely, a harmonic oscillator).
We want briefly to reiterate that this analog fault model is not the state-of-the-art, nor the more sophisticated way to describe the true behavior of an extended, potentially non planar fault embedded in a continuous medium, where heterogeneous properties can come into the game. Nevertheless, it has the great advantage to make it possible to simulate an entire seismic cycle (i.e., several Myr) of a fault with very limited computational resources and in very short computational times [Tse and Rice, 1986;Boatwright and Cocco, 1996;Gomberg et al., 1997;Bizzarri, 2012c].
An intermediate level of complication (in between the two end members of 1-D system and truly 3-D extended fault models) is represented by the model proposed by Burridge and Knopoff [1967], consisting in multiple masses linked together by springs with different elastic constants to model the spatial heterogeneities of the host rocks [Huang et al., 1992], neglected by definition on the case of a single mass-spring model. In a separate study we will explore the role of variable L in a more realistic 3-D setting, where the two modes of shear rupture (mode II and mode III are coupled together).
We conduct a large class of simulations, where we adopt various values of reference scale length, by assuming typical values expected for real-world fault [e.g., Marone and Kilgore, 1993]. In our models we do not change the value of the normal stress n eff , although we know that dynamic perturbations of the effective normal stress acting on the fault (as well as bi-material effects, free surface interactions, thermal pressurization of pore fluids) can significantly affect the energy balance and thus ultimately the rupture dynamics, [e.g., Nielsen et al., 1998;Oglesby et al., 1998;2000;Bizzarri and Cocco, 2006a;2006b;Rubin and Ampuero, 2007;Scala et al., 2017;2019]. Indeed, temporal variations of n eff , potentially coupled with variations of L, can further contribute to affect the stability condition of the system, since both parameters enters in the critical stiffness k c (see equation (10)). Due to the non linear nature of the system, even in the simplified 1-D idealization assumed here, it is not possible to a priori predict how these mutual variations can change the behavior of the system. Incidentally, we recall here that the effect of variation of n eff have been reported in Bizzarri [2012b] again with a mass-spring model.
The numerical experiments presented in this work clearly demonstrate that the introduction of a variable L significantly affects the dynamics of a fault. Namely, we observe that, with a variable scale length: I. the cyclicity of the system is not destroyed, in that the slider is able to enter in its limit cycle and therefore the recurrence time remains constant through time; II. with respect to the reference case (where L does not vary and equals L ref over the whole duration of the seismic cycle) the events have lower magnitude, in that the amount of slip developed and during each instability is lower compared to the reference case; III. the stress release (namely, the difference between the peak traction and the residual traction after the total release) is reduced compared to that typical of the reference configuration. This result is important, since the stress release is connected to the emitted energy. Indeed, we know that in an extended model the size of the cohesive zone (i.e., the region of the fault where the stress drop is accomplished) is directly controlled by the length scale (either d 0 of L in the contexts of the slip-dependent laws of in the framework of the rate and state friction). Moreover, it is also well known that that the increasing slip rate is connected to the rupture at which the rupture front advances and to a shrinking of the cohesive itself [e.g., Bizzarri and Das, 2012] and this finally reflects into an increase of the frequency content of the emitted seismic radiation. At the present knowledge, it is not clear how potential temporal variations of L during the propagation of an extended rupture can affect the size of the cohesive zone and thus the seismic energy. This would be the subject of a separate paper in preparation; IV: the minimum values of the slip velocity attained at the end of the deceleration phase are higher; V: the inter-event time T cycle (see equation (8)) is reduced. The variation can even exceed 50 % of the reference cycle time; this is a natural consequence of previous results iii) and iv); VI: the velocity-weakening regime (b > a) and the initially unstable behavior ( < 1) is not more a sufficient condition to obtain dynamic instabilities; the variations of L can produce > 1 during the motion of the slider and therefore, deepening on L ref and on the assumed governing model, the accelerating phase can be impeded and failure event can abort.
These conclusions hold also for other values of model parameters, such as a and b (not shown here for brevity).
Our numerical simulations show very clearly the behavior of the system. The external load causes a gradual increase of the sliding speed , which in turn affects L, that starts to becomes larger than L ref . This time level is indentified by the vertical lines in Figure 9 (red and blue color for models (3) and (4), respectively). Starting from this time the evolution of the state variable departs from its natural, reference behavior (see Figure 9b). The system then undergoes toward the instability, by exceeding the quasi-static limit of v c = 0.01 mm/s (see Figure 9a); it enters into the dynamic regime, characterized by an accelerating phase. This causes the most relevant variation of L (see Figure 9c). This variation changes the value of the critical stiffness (equation (10)) and can even cause (equation (9)) to be greater than unity. In other words, the accelerating system -which started from an inherently unstable regime, due to the assumed model and constitutive parameters -spontaneously enters into the stable range ( > 1; see section 3.3). At this time v > v CS (equations (11) and (12) for models (3) and (4), respectively); this time instant is denoted by big symbols in the inset of Figure 9c. Due to the high degree of instability of the system (guaranteed by the assumed RD model), in this case the accelerating phase continues, the slip velocity reaches its peak and finally starts to decelerate. This results in a decrease in L (see Figure 9c) and consequently in . At a given time the stiffness ratio becomes again lower than unity (or correspondingly the slip velocity returns again to values lower than v CS ); this time level is denoted by the right vertical segments of the rectangles in Figure 9. Finally, L equals L ref again, because the fast deceleration of the system makes v very low (Figure 9a), and the state variable dumps into the new steady state condition (panel b); see also Bizzarri and Cocco [2003]. The coseismic phase is completed and then the intersesimic phase begins, leading to a further instability event.
In this paper we have seen that the time variable L does not prevent the fault to reproduce the so-called characteristic earthquake (CE). This concept, which of course is inherently attractive, is routinely employed in probabilistic hazard analysis and in building codes. Moreover, although its physical meaning has been scrutinized and widely criticized [Bizzarri, 2012c], the concept of recurrence time often founds the framework of earthquake prediction. Indeed, it has been demonstrated that several physical phenomena can destroy the concept of CE; the pore fluids thermal pressurization [Mitsui and Hirahara, 2009], wear processes [Bizzarri, 2010], porosity and permeability evolution [Mitsui and Cocco, 2010;Bizzarri, 2012b], time variations of parameters a and b due to temperature-dependence [Bizzarri, 2011b]. Although the cyclicity of the system survives with both model (3) and (4), the variations of L shed light on the fact that it would be impossible to a priori predict the exact value of the recurrence time. Indeed, a recent study [Barbot et al., 2012] suggested an analytical expression for T cycle , which is postulated to depend on coseismic (v co ) and interseismic (v inter ) slip velocities: While v co can be estimated roughly in 1 m/s, v inter is unpredictable and strongly depends on the value of L ref and on the assumed friction model, as we have clearly shown.
To summarize, we have demonstrated that, independently of the assumed parameters and on the adopted friction model, the velocity-dependence of the scale length is able to produce significant modifications to the reference behavior even of a simple, idealized fault model, such as the spring-slider system.