WANG_correctedOK:Layout 6 ANNALS OF GEOPHYSICS, 56, 3, 2013, R0332; doi:10.4401/ag-5548 R0332 Stability analysis of slip of a one-body spring-slider model in the presence of thermal pressurization Jeen-Hwa Wang Institute of Earth Sciences, Academia Sinica, Taipei, Taiwan ABSTRACT 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 llcr, thus not resulting in unstable motions. The phase portraits of v versusd for l<30 MPa/m show the existence of a non-zero unstable fixed point. 1. Introduction After an earthquake ruptures, the shear stress, x, on the fault plane decreases from a static xo to a dy- namic xd, and finally becomes xf. The weakening mech- anism 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 ener- getically equivalent, then it would be virtually impossi- ble 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 state- dependent law [Bizzarri and Cocco 2005, Wang 2009a, Wang 2012]. Bizzarri and Cocco [2003] showed that the rate-weakening law can also yield slip-weakening be- havior of shear-stress evolution. Stability analysis of slip of a one-body (one-degree- of-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 sta- bility 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 equilib- rium 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 rate- dependent 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 consid- ered, 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 and Tullis 2002, Di Toro et al. 2004], and melting lubrication [Sib- son 1975, Spray 1993, Tsutsumi and Shimamoto 1997, Hirose and Shimamoto 2005], among others. A sum- mary of the mechanisms can be found in Rice [2006] and Bizzarri [2009]. Except for flash heating and melt- ing lubrication, heat and the pore fluid pressure both have major roles in reducing the frictional stress. Ther- Article history Received January 6, 2012; accepted May 28, 2013. Subject classification: Slip stability, One-body spring-slider system, Friction, Thermal pressurization, AUD and SOP models. mal pressurization has long been considered to be an important mechanism in controlling fault dynamics. Rice [2006] proposed two end-member models for ther- mal pressurization: the adiabatic-undrained-deformation (AUD) model, and the slip-on-a-plane (SOP) model. Wang [2009b] showed that the AUD model is more ap- propriate than the SOP model to interpret the observed xf – d function. To understand slip behavior of a fault under slip- dependent friction caused by thermal pressurization, it is significant to analyze the stability of slip in the neigh- borhood of the equilibrium state on the basis of a one- body 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 adi- abatic 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 xf – d functions will be carried out. For this study, the stability condition for the SOP model will also be de- scribed. Second, the stability of slip at and near v=0 will be analyzed on the basis of a set of equations to repre- sent the motions of the one-body model in the pres- ence of AUD-type thermal pressurization. Third, numerical simulations will be performed from those equations. 2. 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, ct=dc/dt is the shear rate, t is the density of the fault zone, tf is the density of the fluids, c is the heat capacity, at is the thermal diffusivity, ah is the hydraulic diffusivity, and np is the inelastic (or plastic) porosity after deformation. The parameter b is the volumetric pore fluid storage coefficient and it is equal to n(bf+bn), where n is the starting porosity before deformation, bf is the isothermal compressibility of the pore fluid (dtf/tf=bfdp), and bn is the isothermal compressibility of the pore space (dn/n=bfdp). is the undrained pres- surization factor that can be written as (mf−mn)/(bf+bn), where mf (dtf/tf=−mffdT) and n (dn/n=mndT) are the isobaric, volumetric thermal expansion coefficients for pore fluids and pore space, respectively. The parame- ters at and ah are equal to K/tc and k/hfb, respectively, where K, k, and hf are the thermal conductivity, per- meability, 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 vn on a spatial scale of the sheared layer that is broad enough to ef- fectively 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 trans- ferred outwards from the fault plane. Before slipping, the effective static shear stress x on the fault plane can be represented by xo=f(vn−po), where f is the static friction coefficient and vn and po are the normal stress and pore fluid pressure, respectively, on the sliding plane (y=0). The xf – d functions caused by thermal pressuriza- tion [Rice 2006] are: (1) For the AUD model, (3) In Equation (3), dc=tch/fK, where h is the thick- ness of the slip zone. (2) For the SOP model, (4) JEEN-HWA WANG 2 T/ t / 1/ / yy /c c c Tt t2 2 2 2 22vc t t t a- =^ ^h h6 @ / t T/ t / t 1/ / y / y ,p /p n ff h p2 2 2 2 2 2 2 2 22b b b t t a K- + = = ^ ^ ^ h h h6 @ Figure 1. The functions of shear stress, x, in terms of slip, d: curves for xf versus d: and lines for xe versus d. xo (=ldL) denotes the effec- tive shear strength. / .expfAUD o cx d x d d= -^ ^h h /L /L .erfcexpf o / SOP 1 2 x d x d d= ) )^ ^ ^h h h 3 To derive Equation (4), the sliding velocity, v (=dd/dt), is set to be constant. In Equation (4), erfc(z) is the com- plementary error function, and L*=4(tc/K)2(at 1/2 + + ah 1/2)2/f 2v. The two parameters dc and L * control the xf – d functions of individual end-member models. In Equa- tion (3), dc is a characteristic displacement of the xf – d function. For Equation (4), there is no characteristic dis- placement for the xf – d function [Rice 2006]. An ex- ample to show the xf – 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 xf with increasing d for each of these models. 3. 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, vn. 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 equi- librium position, along the horizontal axis. The slider is driven by elastic traction from the spring, which moves with a velocity Vp. The slip at the loading point is do=Vpt. The elastic traction exerted on the slider by the spring is: (5) The slider is also subjected to slip-dependent fric- tional stress, xf (d). Hence, the equation of motion of the system is: (6) 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 ad- dition, the AUD model represents conditions over a suf- ficiently short time scale to neglect fluid and heat dif- fusion during the earthquake rupture, while the SOP model represents those over a longer time scale, to neg- lect 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 ap- proach. Inserting Equation (3) into Equation (6) leads to: (7) 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(d2d/dt2)=0. At this moment, do is dL, thus leading to ldL−xo=0 from Equation (7). This gives dL=xo/l. When do is slightly greater than dL, 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 earth- quake is several orders shorter than the loading time, dL can be considered to be constant during faulting. After the slider moves, xe decreases with increasing d. The relationship of xe=l(do−d) is shown in Figure 1 by the straight line with slope −l. This line crosses the ver- tical axis for x and the horizontal axis for d at x=xo=ldL and d=dL, respectively. 4. Analytical study The xf – d function is shown by the bold curve in Figure 1. The derivatives of xf with respective to d for the AUD model are: (8) and for the SOP model they are: (9) Equations (8) and (9) lead to dxf/dd=−xo/dc and dxf/dd→∞ at d=0, respectively. lcr=|dxf/dd| at d=0 is defined as the characteristic stiffness of the system. Hence, lcr equals xo/dc for the AUD model, and ∞ for the SOP model. For the AUD model, the xe – d func- tion with l>lcr is shown with a line marked by l>lcr in Figure 1a. Obviously, the line cannot cross the xf – d function, and xe is always lcr in Figure 1. This line does THERMAL PRESSURIZATION ON SLIP STABILITY d /dt / .expe coo 2 2 d d dt d l d x= -- -^ ^ ^h h h d /d /L erfc /L / L .exp 1 f o / /1 2 1 2 x d x d d r d = -= ) ) )^ ^ ^h h h6 @ Figure 2. A one-body dynamical spring-slider model: m, mass of the slider; vn, normal stress; Vp, driving velocity; xf, frictional stress; l, spring constant; d, slip of the slider; and do, slip of the loading point. .e ox d l d d= -^ ^h h d /dt d dtv/ ,e f 2 2t d t x x= = -^ ^h h d /d / /expf o c cx d x d d d=- -^ h cross the xf – d function in Figure 1a, although not in Figure 1b. From d=0 to the displacement related to the intersection point, xe is higher than xf, thus indicating unstable motion. Hence, lcr is the critical stiffness (spring constant) for controlling stable or unstable mo- tions of the system. The ratio xo/dc is a function of the thermal and hydromechanical parameters of a fault zone; Wang [2009b] indicated that the two microscopic physical pa- rameters xo and dc can be evaluated from macroscopic seismological observations. When l>lcr=xo/dc, 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 lxo/dc, and weak coupling when l>T2. The normalized equations of motions from Equa- tion (6) are: (10a) (10b) (10c) From Taylor’s series, the first-order approximation of e−x is 1−x as |x|<<1; i.e., e−x�1−x. This approxi- mation makes Equation (10b) become: (11) Differentiating Equation (11) leads to: (12) Inserting Equations (10a) and (10c) into Equation (12) results in: (13) The Laplace transformation of Equation (13) is: (14) where Y(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)=s2−[(xo/dc)− −l]/t=0. The slipping state is stable if the equation D(s)=0 has no solutions so with Re[so]>0, and unstable if the equation D(s)=0 has some solutions so with Re[so]>0. The roots of Equation (15) are: (16) If llcr, 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: llcr represent the unstable (supercritical), marginally stable (critical), JEEN-HWA WANG 4 y V / z e T ;//z e p o x t x 2t x= - -= - -^ ^ ^h h h s Y s T T s T1/ sT 1/ 1/ Y /2 3 1 3 2= + -^ ^ ^ ^h h h h6 @ Y s T s s T T T s s 1/T / 1/ 1/ / / / / / .o c 2 3 2 1 3 2 2l t x d l t = - - = = - - ^ ^ ^ ^ ^ h h h h h 6 6 @ @" , s / / / .o c / cr /1 2 1 2! !x d l t l l t= - = -! ^ h6 6@ @ @ @ x y/ /V y/T ;t c p 1d= =^ h z 1 y y/ / V 1 /T .t o p 3x l= - = -^ ^ ^h h h y z 1 x /T .t 2= - +^ h y z x /T .t t tt 2= +^ h y 1 y y T T/T / / .tt 3 1 2= - +^ h6 @ 5 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 above- mentioned numerous articles [Gu et al. 1984]. 6. Numerical simulations Simulations of slip of the one-body system in the presence of AUD-type thermal pressurization are per- formed on the basis of Equation (10), using the fourth- order Runge-Kutta algorithm [Press et al. 1986]. As mentioned above, f represents the static friction coeffi- cient between the bottom of the slider and the moving plate. When vn−po and f are set to be 100 MPa and 0.6, respectively, the effective shear stress, i.e., xo=f(vn−po), is 60 MPa. The value of dc is 2 m. Hence, the critical stiffness, lcr, 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/m1. In Figure 3c, the bisection line intersects two points. Hence, there are two fixed points: one is the zero point (denoted by P0=(0.0, 0.0)), and the other is a nonzero point (denoted by P1). Obviously, for the two models, the values of h at P0 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 P1 is negative, with a magnitude >1. Hence, the nonzero point is an unstable fixed point. As the velocity reaches its maximum, xe−xf becomes neg- ative, thus resulting in resistance to decelerate the slider. Then the magnitude of xe−xf reduces with de- creasing velocity, and finally becomes zero when the slider stops. In other words, the portrait returns to P0. The driving force continues to pull the slider, thus in- creasing the value of xe on it. When xe > xf, the slider moves again, and thus the portrait in Figure 3c repeats. The phase portraits of v/vmax versus / max are plotted in Figure 3d. Essentially, v/vmax increases and then decreases with increasing d/dmax. The dashed line in Figure 3c represents the bisection line, which inter- sects the phase portrait at two points: one at the zero point (P0) and the other at a nonzero point (P1). The values of h at P0 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 P1, 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 nu- merous cases, with l=5, 10, 15, 20, 25, 29, 30, and 31 MPa/m. The first six cases with llcr, and marginally stable (or critical) when l=lcr, and unstable (or supercritical) when l> /ColorImageDict << /QFactor 0.15 /HSamples [1 1 1 1] /VSamples [1 1 1 1] >> /JPEG2000ColorACSImageDict << /TileWidth 256 /TileHeight 256 /Quality 30 >> /JPEG2000ColorImageDict << /TileWidth 256 /TileHeight 256 /Quality 30 >> /AntiAliasGrayImages false /CropGrayImages true /GrayImageMinResolution 150 /GrayImageMinResolutionPolicy /OK /DownsampleGrayImages true /GrayImageDownsampleType /Bicubic /GrayImageResolution 300 /GrayImageDepth -1 /GrayImageMinDownsampleDepth 2 /GrayImageDownsampleThreshold 1.10000 /EncodeGrayImages true /GrayImageFilter /DCTEncode /AutoFilterGrayImages true /GrayImageAutoFilterStrategy /JPEG /GrayACSImageDict << /QFactor 0.15 /HSamples [1 1 1 1] /VSamples [1 1 1 1] >> /GrayImageDict << /QFactor 0.15 /HSamples [1 1 1 1] /VSamples [1 1 1 1] >> /JPEG2000GrayACSImageDict << /TileWidth 256 /TileHeight 256 /Quality 30 >> /JPEG2000GrayImageDict << /TileWidth 256 /TileHeight 256 /Quality 30 >> /AntiAliasMonoImages false /CropMonoImages true /MonoImageMinResolution 1200 /MonoImageMinResolutionPolicy /OK /DownsampleMonoImages true /MonoImageDownsampleType /Bicubic /MonoImageResolution 1200 /MonoImageDepth -1 /MonoImageDownsampleThreshold 1.08250 /EncodeMonoImages true /MonoImageFilter /CCITTFaxEncode /MonoImageDict << /K -1 >> /AllowPSXObjects false /CheckCompliance [ /None ] /PDFX1aCheck false /PDFX3Check false /PDFXCompliantPDFOnly false /PDFXNoTrimBoxError true /PDFXTrimBoxToMediaBoxOffset [ 0.00000 0.00000 0.00000 0.00000 ] /PDFXSetBleedBoxToMediaBox true /PDFXBleedBoxToTrimBoxOffset [ 0.00000 0.00000 0.00000 0.00000 ] /PDFXOutputIntentProfile (None) /PDFXOutputConditionIdentifier () /PDFXOutputCondition () /PDFXRegistryName (http://www.color.org) /PDFXTrapped /Unknown /CreateJDFFile false /SyntheticBoldness 1.000000 /Description << /ENU (Use these settings to create PDF documents with higher image resolution for high quality pre-press printing. The PDF documents can be opened with Acrobat and Reader 5.0 and later. These settings require font embedding.) /JPN /FRA /DEU /PTB /DAN /NLD /ESP /SUO /NOR /SVE /KOR /CHS /CHT /ITA >> >> setdistillerparams << /HWResolution [2400 2400] /PageSize [595.000 842.000] >> setpagedevice