Stability analysis of slip of a one-body spring-slider model in the presence of thermal pressurization

This study is focused on the stability of slip in a one-body spring-slider system, with stiffness l of the slider, in the presence of slip-dependent friction caused by thermal pressurization for two end-member models of thermal pressurization; i.e., the adiabatic-undrained-deformation (AUD) model and the slip-on-a-plane (SOP) model. Analytical studies based on the functions of frictional stress, xf , versus slip, d, of the two models show that lcr=|dxf/dd| at d=0 is the critical stiffness of the system. lcr is a finite positive value for the AUD model, and infinity for the SOP model. Slip is stable when l>lcr and unstable when l<lcr. The analytical stability conditions of slip from the equations of motions of a onebody spring-slider model with the xf – d function of the ADU model are also those from |dxf/dd| at d=0. Under ADU thermal pressurization, the modeled final slip and maximum velocity for l=5 to 31 MPa/m, with lcr=30 MPa/m, both decrease with increasing l. The motions become very weak when l>lcr, thus not resulting in unstable motions. The phase portraits of v versus d for l<30 MPa/m show the existence of a non-zero unstable fixed point.


Introduction
After an earthquake ruptures, the shear stress, x, on the fault plane decreases from a static x o to a dynamic x d , and finally becomes x f .The weakening mechanism is complex [Bizzarri 2009] and can be either slip-dependent, or rate-dependent and state-dependent [Dieterich 1979, Ruina 1983, Marone 1998, Wang 2002, Bizzarri and Cocco 2006a].Bizzarri [2011] showed that if two faults that obey different friction laws are energetically equivalent, then it would be virtually impossible to distinguish between them, even on the basis of their frequency content and synthetic ground motions.Note that slip, d, and the rate or sliding velocity, v, are both a function of time.For a long time, there has been debate concerning the constitutive laws.Some studies are in favor of the slip-dependent law [Ohnaka 2004], while other studies prefer the rate-dependent and statedependent law [Bizzarri andCocco 2005, Wang 2009a, Wang 2012].Bizzarri and Cocco [2003] showed that the rate-weakening law can also yield slip-weakening behavior of shear-stress evolution.
Stability analysis of slip of a one-body (one-degreeof-freedom) spring-slider model (see Figure 1) is usually carried out to understand the behavior of faulting.For rate-dependent and state-dependent friction, slip stability has been analyzed in numerous studies [Rice 1983, Ruina 1983, Rice and Ruina 1983, Gu et al. 1984, Blanpied and Tullis 1986, Rice and Tse 1986, Gu and Wong 1991, Belardinelli and Belardinelli 1996, He et al. 1998, Reinen 2000, Ryabov and Ito 2001].The equilibrium state of the one-body system subjected to friction exists when the slider moves under a driving velocity, and in the neighborhood of this state it shows that the model parameters control the stability of slip of the system.As mentioned above, stability analyses of this type of system have mainly been focused on the ratedependent and state-dependent friction law.However, this kind of analysis has not yet been performed for the case of slip-dependent friction.
To interpret the decrease in the effective frictional stress, several physical mechanisms have been considered, including: hydrodynamic lubrication [Brodsky and Kanamori 2001], thermal pressurization [Sibson 1973, Fialko 2004, Bizzarri and Cocco 2006b, 2006c, Rice 2006, Wang 2009b], flash faulting [Rice 1999, Tullis and Goldsby 2003a, Tullis and Goldsby 2003b, Prakash 2004, Prakash and Yuan 2004, Hirose and Shimamoto 2005, Rice 2006], gel formation [Goldsby andTullis 2002, Di Toro et al. 2004], and melting lubrication [Sibson 1975, Spray 1993, Tsutsumi and Shimamoto 1997, Hirose and Shimamoto 2005], among others.A summary of the mechanisms can be found in Rice [2006] and Bizzarri [2009].Except for flash heating and melting lubrication, heat and the pore fluid pressure both have major roles in reducing the frictional stress.Ther-mal pressurization has long been considered to be an important mechanism in controlling fault dynamics.Rice [2006] proposed two end-member models for thermal pressurization: the adiabatic-undrained-deformation (AUD) model, and the slip-on-a-plane (SOP) model.Wang [2009b] showed that the AUD model is more appropriate than the SOP model to interpret the observed x f -d function.
To understand slip behavior of a fault under slipdependent friction caused by thermal pressurization, it is significant to analyze the stability of slip in the neighborhood of the equilibrium state on the basis of a onebody spring-slider model subject to slip-dependent shear stress caused by thermal pressurization.Unlike rate-dependent and state-dependent friction, to date there is a lack of this a kind of study.As the SOP model is defined based on a non-zero constant slip velocity, v, this model cannot work at v=0.Hence, in this study, an attempt is made mainly to investigate slip stability of a one-body spring-slider model subject to slip-dependent shear stress caused by thermal pressurization under adiabatic and undrained conditions; i.e., the AUD model.For this purpose, three approaches will be used.First, an analytic study of stability conditions based on the x f -d functions will be carried out.For this study, the stability condition for the SOP model will also be described.Second, the stability of slip at and near v=0 will be analyzed on the basis of a set of equations to represent the motions of the one-body model in the presence of AUD-type thermal pressurization.Third, numerical simulations will be performed from those equations.

One-dimensional thermal pressurization model
Considering a one-dimensional fault model, the x and y axes denote the directions along and normal to the fault, respectively.Under thermal pressurization, the energy and fluid mass conservation equations can be written as [Rice 2006]: (1) (2) where t is the time, T is the temperature, c is the shear strain, c t =dc/dt is the shear rate, t is the density of the fault zone, t f is the density of the fluids, c is the heat capacity, a t is the thermal diffusivity, a h is the hydraulic diffusivity, and n p is the inelastic (or plastic) porosity after deformation.The parameter b is the volumetric pore fluid storage coefficient and it is equal to n(b f +b n ), where n is the starting porosity before deformation, b f is the isothermal compressibility of the pore fluid (dt f /t f =b f dp), and b n is the isothermal compressibility of the pore space (dn/n=b f dp). is the undrained pressurization factor that can be written as (m f −m n )/(b f +b n ), where m f (dt f /t f =−m f fdT) and n (dn/n=m n dT) are the isobaric, volumetric thermal expansion coefficients for pore fluids and pore space, respectively.The parameters a t and a h are equal to K/tc and k/h f b, respectively, where K, k, and h f are the thermal conductivity, permeability, and fluid viscosity, respectively.Rice [2006] proposed two end-member models for thermal pressurization: the AUD model and the SOP model, detailed explanations of which can be found in Rice [2006].Only a brief description is given here.The AUD model corresponds to a homogeneous simple shear strain c at a constant normal stress v n on a spatial scale of the sheared layer that is broad enough to effectively preclude heat or fluid transfer.The SOP model shows that all sliding with a constant velocity is merely on the fault plane and heat is simultaneously transferred outwards from the fault plane.Before slipping, the effective static shear stress x on the fault plane can be represented by x o =f(v n −p o ), where f is the static friction coefficient and v n and p o are the normal stress and pore fluid pressure, respectively, on the sliding plane (y=0).
The x f -d functions caused by thermal pressurization [Rice 2006] are: (1) For the AUD model, In Equation (3), d c =tch/fK, where h is the thickness of the slip zone.

/ . exp
To derive Equation ( 4), the sliding velocity, v (=dd/dt), is set to be constant.In Equation (4), erfc(z) is the complementary error function, and The two parameters d c and L * control the x f -d functions of individual end-member models.In Equation (3), d c is a characteristic displacement of the x f -d function.For Equation ( 4), there is no characteristic displacement for the x f -d function [Rice 2006].An example to show the x f -d function is given in Figure 1a for the AUD model, and Figure 1b for the SOP model.Figure 1 thus shows the decrease in x f with increasing d for each of these models.

One-body spring-slider model
The dynamic one-body spring-slider model (see Figure 2) consists of a slider of mass m that is linked by a spring with stiffness l.The slider is also exerted by a normal stress, v n .The bottom area of the slider is A. At time t=0, the slider rests in an equilibrium state.Its location is d, measured with respect to its initial equilibrium position, along the horizontal axis.The slider is driven by elastic traction from the spring, which moves with a velocity V p .The slip at the loading point is d o =V p t.The elastic traction exerted on the slider by the spring is: (5) The slider is also subjected to slip-dependent frictional stress, x f (d).Hence, the equation of motion of the system is: where t is the mass per unit length, with a unit of kg/m.
As mentioned above, the SOP model is defined on the basis on a nonzero constant slip velocity, v, and thus this model cannot be used for dynamic analysis.In addition, the AUD model represents conditions over a suf-ficiently short time scale to neglect fluid and heat diffusion during the earthquake rupture, while the SOP model represents those over a longer time scale, to neglect shear zone thickness.Hence, the stability analysis of the AUD model is valid for seismic events with short slip durations.On the other hand, the SOP model cannot work for such events.Hence for the AUD model, only Equation ( 3) is taken into account for the dynamic approach.Inserting Equation (3) into Equation ( 6) leads to: Immediately before the slider starts to move, the displacement is d=0, and the total stress exerted on the slider is zero; i.e., t(d 2 d/dt 2 )=0.At this moment, d o is d L , thus leading to ld L −x o =0 from Equation ( 7).This gives d L =x o /l.When d o is slightly greater than d L , the elastic traction from the extended spring can overcome the effective breaking strength, and then it pushes the slider to move.As the rupture duration of an earthquake is several orders shorter than the loading time, d L can be considered to be constant during faulting.After the slider moves, x e decreases with increasing d.The relationship of x e =l(d o −d) is shown in Figure 1 by the straight line with slope −l.This line crosses the vertical axis for x and the horizontal axis for d at x=x o =ld L and d=d L , respectively.

Analytical study
The Equations ( 8) and ( 9) lead to dx f /dd=−x o /d c and dx f /dd→∞ at d=0, respectively.l cr =|dx f /dd| at d=0 is defined as the characteristic stiffness of the system.Hence, l cr equals x o /d c for the AUD model, and ∞ for the SOP model.For the AUD model, the x e -d function with l>l cr is shown with a line marked by l>l cr in Figure 1a.Obviously, the line cannot cross the x f -d function, and x e is always <x f , thus resulting in stable motion.This situation does not exist for the SOP model, because l cannot be larger than l cr =∞.When is smaller than l cr =x o /d c for the AUD model or l cr =∞ for the SOP model, the x e -d function is shown with the line marked by l>l cr in Figure 1.This line does  cross the x f -d function in Figure 1a, although not in Figure 1b.From d=0 to the displacement related to the intersection point, x e is higher than x f , thus indicating unstable motion.Hence, l cr is the critical stiffness (spring constant) for controlling stable or unstable motions of the system.
The ratio x o /d c is a function of the thermal and hydromechanical parameters of a fault zone; Wang [2009b] indicated that the two microscopic physical parameters x o and d c can be evaluated from macroscopic seismological observations.When l>l cr =x o /d c , the stiffness between the driving force and the slider is higher than the characteristic stiffness of the system, and thus the system shows stable motions.On the other hand, when l<l cr =x o /d c , the stiffness between the driving force and the slider is lower than the characteristic stiffness of the system.This leads to unstable motions of the system.The driving force on a fault system comes mainly from the moving plate.Hence, x o /d c is an important factor in controlling the degree of coupling between a moving plate and a fault system: strong coupling when l>x o /d c , and weak coupling when l<x o /d c .The results suggest that large earthquakes can occur more easily in areas with weaker coupling than those with stronger coupling.
At the intersection point of the x e -d and x f -d functions, x e equals x f , thus leading to d 2 d/dt 2 =0.In Figure 1, the values of d c for the AUD model and L * for the SOP model are equal.However, under the same x o , the value of d at the intersection point is larger for the AUD model than for the SOP model.Although the total stress is zero at the intersection point, the slider can still move forwards because the velocity of the slider is greater than zero.This could result in overshoot of the system and will be discussed below.

Stability analysis
To study the stability of slip of the one-body system in the presence of AUD-type thermal pressurization under adiabatic and undrained conditions, Equation ( 6) is normalized by defining x=d/d c , y=v/V p , and z=x e /x o .The temporal derivatives of x, y and z are x t =dx/dt=(dd/dt)/d c =v/d c , y t =dy/dt=(dv/dt)/V p = =(d 2 d/dt 2 )/V p , and z t =dz/dt= (dx e /dt)/x o =[d(V p t− −d)/dt]/x o =[lV p −l(dd/dt)]/x o , respectively.In addition, T 1 =d c /V p , T 2 =tV p /x o , and T 3 =x o /lV p represent three generalized periods of the system.In general, the values of the related parameters are: l=~10 MPa/m and t=~10 7 kg/m [Gu and Wong 1991] [Rice 2006], f=~0.6, and V p =~10 −9 m/s.This leads to T 1 =~10 9 s, T 2 = =~10 −3 s, and T 3 =~10 9 s.Obviously, TT 3 >>T 2 .
The normalized equations of motions from Equa-tion (6) are: From Taylor's series, the first-order approximation of e −x is 1−x as |x|<<1; i.e., e −x 1−x.This approximation makes Equation (10b) become: If l<l cr , s is a real number.The positive root leads to unstable motions, while the negative root results in stable motions.If l=l cr , s is null.This root represents the stable inflected node in the phase portrait, and slip is marginally stable.Hence, the roots represent the saddle point in the phase portrait.If l>l cr , s is an imaginary number.This always results in stable oscillations.Hence, the positive and negative imaginary roots denote the center of oscillations in the phase portrait.Consequently, the three conditions: l<l cr , l=l cr , and l>l cr represent the unstable (supercritical), marginally stable (critical), and stable (subcritical) states of the system, respectively.This is consistent with the analytical results shown in the last section.Meanwhile, the three conditions are in agreement with the counterparts in the rate-dependent and state-dependent friction law, as shown in the abovementioned numerous articles [Gu et al. 1984].

Numerical simulations
Simulations of slip of the one-body system in the presence of AUD-type thermal pressurization are performed on the basis of Equation ( 10), using the fourthorder Runge-Kutta algorithm [Press et al. 1986].As mentioned above, f represents the static friction coefficient between the bottom of the slider and the moving plate.When v n −p o and f are set to be 100 MPa and 0.6, respectively, the effective shear stress, i.e., x o =f(v n −p o ), is 60 MPa.The value of d c is 2 m.Hence, the critical stiffness, l cr , is 30 MPa/m.The density is t=10 7 kg/m.First, the slip-dependent variation in x is simulated for an unstable case, with l=10 MPa/m<l cr .The plots of x versus d are shown in Figure 3a Here, v reaches the maximum value and thus the slider can still move forward.Beyond this point, the frictional stress is higher than the elastic traction, and thus the slider slows down.Figure 3a shows that the difference between x e and x f , i.e., x e −x f , becomes negative, and thus forming resistance to decelerate the slider.the slider slow down and finally stop when d=d F .
The temporal variations in normalized acceleration (a/a max ), normalized velocity (v/v max ), and normalized displacement (d/d max ) are plotted in Figure 3b.The acceleration and velocity first increase to their respective peaks, and then decrease.The acceleration is zero when the velocity reaches its peak and becomes negative due to a negative value of x e −x f , as mentioned above.This makes the slider slow down, and finally the slider stops when the velocity is zero.The displacement increases monotonously from 0 to its peak.
The phase portrait, which is a plot of a physical quantity versus another physical quantity of an object in a dynamical system, can represent the dynamical behavior of the system [Thompson and Stewart 1986].The phase portraits of x/x max (x=x e −x f ) versus v/v max are plotted in Figure 3c.When v/v max increases from 0 to the maximum value, x/x max is positive.Its magnitude increases and then decreases with increasing v/v max .The dashed line in Figure 3c is the bisection line.The intersection point of the bisection line with a phase portrait is called a fixed point [Thompson and Stewart 1986].Let h be the magnitude of slope of the THERMAL PRESSURIZATION ON SLIP STABILITY tangential line at a fixed point.The stability of a fixed point is dependent upon h: it is stable when h<1, marginally stable when h=1, and unstable when h>1.In Figure 3c, the bisection line intersects two points.Hence, there are two fixed points: one is the zero point (denoted by P 0 =(0.0, 0.0)), and the other is a nonzero point (denoted by P 1 ).Obviously, for the two models, the values of h at P 0 are >1 and thus the zero point is a stable fixed point.This is the same as the conclusion from the previous two sections.On the other hand, the value of h at P 1 is negative, with a magnitude >1.Hence, the nonzero point is an unstable fixed point.As the velocity reaches its maximum, x e −x f becomes negative, thus resulting in resistance to decelerate the slider.Then the magnitude of x e −x f reduces with decreasing velocity, and finally becomes zero when the slider stops.In other words, the portrait returns to P 0 .The driving force continues to pull the slider, thus increasing the value of x e on it.When x e > x f , the slider moves again, and thus the portrait in Figure 3c repeats.
The phase portraits of v/v max versus / max are plotted in Figure 3d.Essentially, v/vmax increases and then decreases with increasing d/d max .The dashed line in Figure 3c represents the bisection line, which intersects the phase portrait at two points: one at the zero point (P 0 ) and the other at a nonzero point (P 1 ).The values of h at P 0 are >1.Hence, the zero point is an unstable fixed point.This is the same as the conclusion gained in the previous two sections.At P 1 , h is negative, but its magnitude is >1.Hence, this non-zero point is an unstable fixed point.
To show the effects of a change of l on slip of the system, numerical simulations are conducted for numerous cases, with l=5, 10, 15, 20, 25, 29, 30, and 31  MPa/m.The first six cases with l<l cr are for unstable slip, the seventh case for marginally stable slip, and the last case for stable slip.Figure 4 shows the plots of maximum (final) slip, d max , and maximum velocity, v max , versus l.As the maximum acceleration, a max , varies between positive and negative values, it is not plotted in Figure 4.It is obvious that d max and v max both decrease with increasing l, and the changing rate is higher for d max than for v max , especially at small l.Obviously, the degree of unstable slip increases with decreasing l.The values of d max and v max are very small when l≥l cr .This is because the elastic traction is smaller than the frictional stress, as shown in Figure 1, thus resisting slip.
Numerical results for l=29, 30, and 31 MPa/m are displayed to examine the detailed variations of motions from unstable to stable states when the values of l are close to l cr .Results for the x e −d and x f −d functions, the temporal variations in a/a max , v/v max , and d/d max , and the phase portraits of v/v max versus d/d max are shown in Figure 5a for l=29 MPa/m, Figure 5b for l=30 MPa/m, and Figure 5c for l=31 MPa/m.The values of a max , v max , and d max used for normalization for these three cases are actually those from the case with l=29 MPa/m.For the three cases, the values of d 1 and d F are close to zero, and thus they are not displayed in Figure 5a1, b1, and c1. Figure 5a2 shows the temporal variations in v/v max and d/d max , which are similar to those displayed in Figure 3b for l=10 MPa/m.Nevertheless, unlike Figure 3b, a/a max in Figure 5a2 is negative when v/v max decreases from the maximum value to zero.Numerical tests show that when l is increased, a max changes from positive to negative and the magnitude of a max decreases.Let a maxn be the largest magnitude of negative accelerations and a maxp be that of positive accelerations.The ratio of a maxn to a maxp increases with l, even though a maxn and a maxp both decrease with l.This can be seen from a comparison between Figure 3b and Figure 5a2.Large values of negative accelerations for large l will resist the motion of the slider.This makes d max and v max decrease with increasing l. Figure 5b2 and c2 show that the values of a/a max , v/v max , and d/d max are very small when l≥l cr =30 MPa/m and the duration time is shorter for l=31 MPa/m than for l=30 MPa/m.Like Figure 3c, Figure 5a3 shows a regular, complete cycle of motions.On the other hand, there are only very small values of v/v max and d/d max in Figure 5b3 and c3.Results again confirm that l is the key parameter in controlling slip stability of the system, and the condition of l≥l cr cannot produce unstable motions.
Although there have been numerous theoretical studies of modeling and simulation to date, the applications of thermal pressurization models to real faults are rare, especially for the correlation of slip stability to faulting.Wang [2009b] applied the thermal pressurization models to explain the slip-dependent stress func- tion on the fault plane during an earthquake, and also to study the slip-dependent radiation efficiency, R, in terms of the parameters of the thermal pressurization model.His results can provide us with something significant and are briefly described here.The radiation efficiency, R, is an important parameter that shows the source property.It is strongly affected by the variation in shear stress with slip.

Conclusions
Analytic studies from the functions of frictional stress versus slip show that the value of |dx e /dd| at d=0 represents the critical stiffness (spring constant), l cr , of the system.l cr is f(v n −p o )/d c for the AUD model and ∞ for the SOP model.For the AUD model, slip of the system is stable (or subcritical) when l>l cr , and marginally stable (or critical) when l=l cr , and unstable (or supercritical) when l<l cr ; for the SOP model only unstable slip exists because is always smaller than l cr =∞.Analytic analysis of slip stability of a one-body spring-slider model with the x f −d function caused by ADU thermal pressurization was carried out.The stability conditions are the same as those from |dx f /dd| at d=0.Numerical simulations of slip were made for the one-body system in the presence of AUD thermal pressurization.In the phase portrait of x/x max versus v/v max for l=10 MPa/m, there are two fixed points: one at the zero point and the other at a nonzero point.The former is a stable fixed point, while the latter is an unstable fixed point.In the phase portrait of v/v max ver-sus d/d max for l=10 MPa/m, there are also two fixed points: one at the zero point and the other at a nonzero point.They are both an unstable fixed point.Numerical simulations for l=5 to 31 MPa/m (with l cr =30 MPa/m) show that the final slip and maximum velocity both decrease with increasing l and the motion of the slider is rapidly arrested after it starts to slide when l≥l cr .The phase portrait of v versus d for l=29 MPa/m show that the zero point is a stable fixed point and a nonzero point is an unstable fixed point.The phase portraits of v versus d for l=30 and 31 MPa/m show very small values of v and d, thus showing marginally stable and stable behavior.

Figure 1 .
Figure 1.The functions of shear stress, x, in terms of slip, d: curves for x f versus d: and lines for x e versus d. x o (=ld L ) denotes the effective shear strength.
x f -d function is shown by the bold curve in Figure 1.The derivatives of x f with respective to d for the AUD model are: (8) and for the SOP model they are: (9)

Figure 2 .
Figure 2. A one-body dynamical spring-slider model: m, mass of the slider; v n , normal stress; V p , driving velocity; x f , frictional stress; l, spring constant; d, slip of the slider; and d o , slip of the loading point.
s) and 1/s are the Laplace transforms of y(t) and 1, respectively.Equation (14) gives: (15) As s=0 is the trivial solution, the stability of slip depends only on the solutions of D(s)=s 2 −[(x o /d c )− −l]/t=0.The slipping state is stable if the equation D(s)=0 has no solutions s o with Re[s o ]>0, and unstable if the equation D(s)=0 has some solutions s o with Re[s o ]>0.The roots of Equation (15) are:(16) : a straight line for x e and a curve for x f .At the intersection point of the x e −d and x f −d functions, x e equals x f , thus yielding d 2 d/dt 2 =0.The slip at this point is denoted by d 1 in Figure 3a.

Figure 3 .
Figure 3. (a) The stress-slip function for the AUD model.(b) The temporal variations in a/a max (dashed-dotted line), v/v max (dashed line), and d/d max (solid line), and the horizontal dotted line denotes d 1 /d max .(c) The phase portrait for x/x max versus v/v max and two fixed points: P 0 =(0,0) and P 1 =(0.849,0.846).The slope at P 1 is -1.78.(d) The phase portrait for v/v max versus d/d max and two fixed points: P 0 =(0,0) and P 1 =(0.783,0.787).The slope at P 1 is -1.81.

Figure 4 .
Figure 4. Plots of d max (solid line, open squares) and v max (dashed line, open circles) versus l.The vertical dotted line displays the critical state with l =30 MPa/m.

Figure 5 .
Figure 5. Numerical results to examine the detailed variations of motions from unstable to stable states when the values of l are close to l cr .(a) For l=29 MPa/m.(b) For l=30 MPa/m.(c) For l=31 MPa/m. 1, stress-slip functions; 2, temporal variations in a/a max (dashed-dotted line), v/v max (dashed line), and d/d max (solid line); 3, phase portraits for v/v max versus d/d max .