Multistable Slip of a One-degree-of-freedom Spring-slider 1 Model in the Presence of Thermal-pressurized 2 Slip-weakening Friction and Viscosity 3 4

This study is focused on multistable slip of earthquakes based on a one-degree-of-freedom spring-slider model in the presence of thermal-pressurized slip-weakening friction and viscosity by using the normalized equation of motion of the model. The major model parameters are the normalized characteristic displacement, Uc, of the friction law and the normalized viscosity coefficient, η, between the slider and background plate. Analytic results at small slip suggest that there is a solution regime for η and γ (= 1/Uc) to make the slider slip steadily. Numerical simulations exhibit that the time variation in normalized velocity, V/Vmax (Vmax is the maximum velocity), obviously depends on Uc and η. The effect on the amplitude is stronger due to η than due to Uc. In the phase portrait of V/Vmax versus the normalized displacement, U/Umax (Umax is the maximum displacement), there are two fixed points. The one at large V/Vmax and large U/Umax is not an attractor, while that at small V/Vmax and small U/Umax can be an attractor for some values of η and Uc. When Uc<0.55, unstable slip does not exist. When Uc ≥ 0.55, Uc and η divide the solution domain into three regimes: stable, intermittent, and unstable (or chaotic) regimes. For a certain Uc, the three regimes are controlled by a lower bound, ηl, and an upper bound, ηu, of η. The values of ηl, ηu, and ηu− ηl all decrease with increasing Uc, thus suggesting that it is easier to yield unstable slip for larger Uc than for smaller Uc or for larger η than for smaller η. When Uc<1, the Fourier spectra calculated from simulation velocity waveforms exhibit several peaks, thus suggesting the existence of nonlinear behavior of the system. When Uc>1, the related Fourier spectra show only one peak, thus suggesting linear behavior of the system.


Introduction
The earthquake ruptures consist of three steps: nucleation, dynamical propagation, and arrest.Due to the lack of a comprehensive model, a set of equations to completely describe fault dynamics has not yet been established, because earthquake ruptures are very complicated.Nevertheless, some models, for instance the crack model and dynamical lattice model, have been developed to approach fault dynamics.Several factors will control earthquake ruptures (see Wang, 2016b; and cited references herein), including at least brittle-ductile fracture rheology, normal stress, re-distribution of stresses after fracture, fault geometry, friction, seismic coupling, pore fluid pressure, elastohydromechanic lubrication, thermal effect, thermal pressurization, and metamorphic dehydration.A general review can be seen in Bizzarri (2009).Among the factors, friction and viscosity are two important ones in controlling faulting.Burridge and Knopoff (1967) proposed a one-dimensional spring-slider model (abbreviated as the 1-D BK model henceforth) to approach fault dynamics.Wang (2000Wang ( , 2012) ) extended this model to a two-dimensional version.The two models and their modified versions have been long and widely applied to simulate the occurrences of earthquakes (see Wang, 2008Wang, , 2012; and cited references therein).In the followings, the one-, two-, three-, few-, and many-body models are used to represent the one-, two-, three-, few-, and many-degree-of-freedom spring-slider models, respectively.
The few-body models have been long and widely used to approach faults (Turcotte, 1992) Since the commonly-used friction laws are nonlinear, the dynamical model itself could behave nonlinearly.A nonlinear dynamical system can exhibit chaotic behaviour under some conditions (Thompson and Stewart, 1986;Turcotte, 1992).This means that the system is highly sensitive to initial conditions (SIC) and thus a small difference in initial conditions, including those caused by rounding errors in numerical computation, yields widely diverging outcomes.This indicates that long-term prediction is impossible in general, even though the system is deterministic, meaning that its future behavior is fully determined by their initial conditions, without random elements.This behavior is known as (deterministic) chaos (Lorenz, 1963).
The frictional effect on earthquake ruptures has been widely studied as mentioned above.However, the studies of viscous effect on earthquake ruptures are rare.The viscous effect mentioned in Rice et al. (2001) was just an implicit factor which is included in the evolution effect of friction law.In this work, I will investigate the effects of thermal pressurized slip-weakening friction and viscosity on earthquake ruptures and the generation of unstable (or chaotic) slip based on a one-body model.

One-body Model
Fig. 1 shows the one-body model whose equation of motion is: where m is the mass of the slider, u and v (=du/dt) are, respectively, the displacement and velocity of the slider, u o is the equilibrium location of the slider, K is the spring constant, F is the frictional force between the slider and the background and a function of u or v, and Φ is the viscous force between the slider and the background and a function of v.The slider is pulled by a driving force F D due to the moving plate with a constant driving velocity, v p , through a leaf spring of strength, K. Hence, the driving force is F D =Kv p t and thus u o =v p t.When F D is slightly larger than the static frictional force, F o , friction changes from static friction strength to dynamic one and thus the slider moves.(1942) first emphasized the importance of viscosity on faulting.

Jeffreys
Frictional melts in faults depend on temperature, pressure, water content, and etc. (Turcotte and Schubert, 1982) and can yield viscosity on the fault plane (Byerlee, 1968).Rice et al. (2001) discussed that rate-and state-dependent friction in thermally activated processes allows creep slippage at asperity contacts on the fault plane.Scholz (1990) suggested that the friction melts would present significant viscous resistance to shear and thus inhibit continued slip.However, Spray (1993Spray ( , 1995Spray ( , 2005) ) stressed that the frictional melts possessing low viscosity could generate a sufficient melt volume to reduce the effective normal stress and thus act as fault lubricants during co-seismic slip.His results show that viscosity remarkably decreases with increasing temperature.For example, Wang (2011) assumed that quartz plasticity could be formed in the fault zone when T>300 o C after faulting and it would lubricate the fault plane at higher T and yield viscous stresses to resist slip at lower T. From numerical simulations, Wang (2007Wang ( , 2016bWang ( , 2017) ) stressed the viscous effect on faulting.Noted that several researchers (Knopoff et al., 1973;Cohen, 1979;Xu and Knopoff, 1994;Knopoff and Ni, 2001;Dragoni and Santini, 2015) took viscosity as a factor in causing seismic radiation to reduce energy during faulting.
The viscosity coefficient, υ, of rocks is mainly controlled by temperature, T.An increase in T will yield partial melting of rocks and thus the viscosity coefficient, υ, first is increased, then reaches the largest value at a particular T, and finally decreases with increasing T The relation between υ and T can be described by the following equation (e.g., Turcotte and Schubert, 1982) The description about the physical models of viscosity can be found in several articles (Jaeger and Cook, 1977;Cohen, 1979;Hudson, 1980;Wang, 2016b).A brief description is given below.For many deformed materials, there are elastic and viscous components.The viscous component can be modeled as a dashpot such that the stress-strain rate relationship is: σ=υ(dε/dt) where σ and ε are the stress and the strain, respectively.Two simple models (shown in Fig. 2 Under a constant tensile stress, the strain will increase, without a upper limit, with time for the Maxwell model; while the strain will increases, with a upper limit, with time for the Kelvin-Voigt model.Wang (2016b) assumed that the latter is more appropriate than the former to be applied to the seismological problems as suggested by Hudson (1980).Hence, the Kelvin-Voigt model is taken in this study.To simplify the problem, only a constant viscosity coefficient is considered in a numerical simulation as given below.The viscous stress at the slider is represented by -υv.
However, it is not easy to directly implement viscosity in a dynamical system as used in this study.Wang (2016b) represented the viscosity coefficient in an alternative way.Viscosity leads to the damping of oscillations of a body in viscous fluids.The damping coefficient, η, depends on the viscosity coefficient, υ, and the linear dimension, R, of the body in a viscous fluid.According to Stokes' law, the η of a sphere of radius R in a viscous fluid of υ is η=6πRυ (cf.Kittel et al., 1968).In order to simplify the problem, the damping coefficient is taken in this study.Hence, the viscous force is Φ=ηv.Noted that the unit of η is N(m/s) -1 .

Friction caused by thermal pressurization
Numerous factors can affect friction (see Wang, 2009Wang, , 2016b; and cited references herein).When fluids are present and temperature changes in faults, thermal pressurization will yield resistance on the fault plane and thus play a significant role on earthquake rupture (Sibson, 1973;Lachenbruch, 1980;Chester and Higgs, 1992;Fialko, 2004;Fialko and Khzan, 2005;Bizzari and Cocco, 2006a,b;Rice, 2006;Wang, 2000Wang, , 2006Wang, , 2009Wang, , 2011Wang, , 2013Wang, , 2016b, 017;, 017;Bizzarri, 2010;Bizzarri, 2011a,b).Rice (2006) proposed two end-members models for thermal pressurization: the adiabatic-undrained-deformation (AUD) model and slip-on-a-plane (SOP) model.He also obtained the shear stress-slip functions caused by the two models.The first model corresponds to a homogeneous simple shear strain ε at a constant normal stress σ n on a spatial scale of the sheared layer that is broad enough to effectively preclude heat or fluid transfer.The second model shows that all sliding is on the plane with τ(0)= f(σ n -p o ) where p o is the pore fluid pressure on the sliding plane (y=0).For this second model, heat is transferred outwards from the fault plane.Although the stress τ sop (u) also shows slip-weakening (Wang, 2009), the SOP model is not appropriate in this study because of the request of a constant velocity for this model.
The shear stress-slip functions, τ(u), caused by the AUD model is: The parameters u c is the characteristic displacements associated with the thickness and some physical properties of fault zone.The stress τ aud (u) displays exponentially with u and thus exhibits slip-weakening friction.Based on the AUD model, Wang (2009) proposed a simplified slip-weakening friction law (denoted by the TP law hereafter): F(u)=F o exp(-u/u c ), where F o is the static frictional force, to study seismic efficiency.Wang (2016bWang ( , 2017) ) applied the law to simulate slip of one-body and two-body spring-slider models.Fig. 3 exhibits F(u) versus u for five values of u c , i.e., 0.1, 0.3, 0.5, 0.7, and 0.9 m.The friction force decreases with increasing u and it decreases faster for smaller u c than for larger u c .Meanwhile, the force drop decreases with increasing u c .For small u, exp(-u/u c ) can be approximated by 1-u/u c (Wang, 2016a(Wang, ,b, 2017)).The parameter u c -1 is almost the decreasing rate, γ, of friction force with slip at small u.Small (large) u c is related to large (small) γ.

Predominant Frequency and Period of the System
To conduct marginal analyses of slip of one-body model with friction, Wang (2016b) used the friction law: F(u)=F o -γu.His results show that the natural periods are T o =2π/(K/m) 1/2 when friction and viscosity are excluded and when friction and viscosity are included.Clearly, T n is longer than Τ o .Eq. ( 4) shows that T n increases with η and γ, thus indicating that friction and viscosity both lengthen the natural period of the system.

Normalization of Equation of Motion
Substituting the TP law and the linear viscous law into Eq.(1) leads to To simplify numerical computations, Eq. ( 5) is normalized based on the following When F D =v p t or Γ D =V p τ, Eq. ( 6) is transformed to a set of three first-order differential equations by defining x=U/U c , y=V/V p , and z=-U+V p τ-ηV p y τ (y t =dy/dτ): As x<<1, e -x ≈1-x and thus Eq. ( 7b) can be approximated by y τ ≈(z-1+x)/V p .The condition of x<<1 shows U/U c <<1. Differential of this equation leads to y ττ ≈ (z τ +x τ )/V p , where y ττ =d 2 y/dτ 2 .Substituting Eqs. ( 7a) and (7b) into this equation gives The homogeneous equation of Eq. ( 8) is Let the general solution be y～e λτ .This leads to [λ 2 +ηλ+(1-/U c )]y=0 or The solutions of Eq. ( 10) are The term -η/2 of Eq. ( 11) leads to e -λ/2 which yields attenuation of y.Define D(η,1/U c ) to be η 2 -4(1-1/U c ).As mentioned above, U c -1 is the normalized decreasing rate of friction, Γ, at U=0. Fig. 4 shows the plot of η versus 1/U c and thus exhibits the root structure of the system.Because η>0 and U c >0, only the plot in the first quadrant is present in Fig. 4. The solid line displays the function: Along the line, we have η 2 =4(1-1/U c ), and thus λ±=-η/2.In other word, the roots are equal and real, and thus the solution is a stable inflected node displayed by a solid circle in Fig. 4. As D(η,1/U c )>0 or η 2 >4(1-1/U c ), the roots are both real and negative.
The solution shows no oscillation and thus is a stable node shown by a solid square in Fig. 4. As D(η,1/U c )<0 or η 2 <4(1-1/U c ), the roots are complex with negative real part.
This results in oscillations of exponentially decaying amplitude.The solution is a stable spiral or a stable focus shown by an open circle in Fig. 4.

Numerical Simulations
Let y 1 =U and thus y 2 =dU/dτ.Eq. ( 6) can be re-written as two first-order differential equations: Eq. ( 12) will be numerically solved using the fourth-order Runge-Kutta method (Press et al., 1986).To simplify the following computations, the value of Γ D is set to be a small constant of 10 -5 , which can continuously enforce the slider to move.
A phase portrait, denoted by y=f(x), is a plot of a physical quantity versus another of an object in a dynamical system (Thompson and Stewart, 1986).The intersection point of the bisection line, i.e., y=x, and f(x) is called the fixed point, that is, f(x)=x.If the function f(x) is continuously differentiable in an open domain near a fixed point x f and |f'(x f )|<1, attraction is generated.In other words, an attractive fixed point is a fixed point x f of a function f(x) such that for any value of x in the domain that is close enough to x f , the iterated function sequences, i.e., x, f(x), f 2 (x), f 3 (x),…, converges to x f .An attractive fixed point is a special case of a wider mathematical concept of attractors.Chaos can be generated at some attractors.The details can be seen in Thompson and Stewart (1986) or other nonlinear literature.In the following plots, the intersection points of the bisection line (denoted by a thin solid line) with the phase portrait of V/V max versus U/U max are the fixed points.To explore nonlinear behavior of a system, the Fourier spectrum F[V(Ω k )], where Ω k =k/δτ is the dimensionless angular frequency at k=0, ..., N-1, is calculated for the simulation velocity waveform through the fast Fourier transform (Press et al., 1986).The bifurcation from a predominant period to others will be seen in the Fourier spectra.
Numerical simulations for the time variation in V/V max , the phase portrait of and (d) for η=0.90 when U c =0.20.Obviously, the time variation in V/V max exhibits cyclic oscillations with a particular period when η<η l =0.86 and has intermittent slip with shorter periods when η>η l .Such a phenomenon holds also for η<5.5.are given in Table 1.The figure exhibits a stable regime when η≤η l , an intermittency or transition regime when η l <η≤η u , and unstable regime when η>η u .

Discussion
As mentioned above, the natural period of the one-body system at low displacements is T o =2π/ω o =2π(m/K) 1/2 in the absence of friction and viscosity and for U c =0.7; and (d) for U c =0.9 when η=0.In the first panels, the time variation in V/V max exhibits cyclic behavior and the amplitude of V/V max decreases and the predominant period of signal increases with increasing U c .This is consistent with Eq.
(5) in which T n increases with U c .Although the four phase portraits are almost similar, their size decreases with increasing U c .The second panels exhibit two fixed points: one at V=0 and U=0 and the second one at larger V and larger V.The slope values at the first fixed points decrease with increasing U c , thus suggesting that the fixed point is not an attractor for small U c and can be an attractor for larger U c .The slope values at the fixed points for smaller U c are greater than 1 and thus they cannot be an attractor.The third panel for each case displays the Fourier spectrum.Fourier spectra show that, in addition to the peak related to the predominant frequency, there are numerous peaks associated with higher frequencies.This shows nonlinear behavior caused by nonlinear friction.The frequency related to the first peak decreases with increasing U c .The amplitude of a peak decreases with increasing U c .The amplitude of a peak decreases with increasing Ω for small U c ; while it first increases up to the maximum and then decreases with increasing Ω for large U c .The amplitude of a peak becomes very small when Ω>0.25.(c) for U c =1.15; and (d) for U c =2.0 when η=0.In the first panels, the time variation in V/V max exhibits cyclic behavior and the amplitude of V/V max remarkably decreases with increasing U c when U c >1.In the second panels, the size of phase portrait decreases with increasing U c and there are two fixed points: the first one at V=0 and U=0 and the second one at larger V and larger V.With comparison to the phase portrait of U c =1.0, the phase portrait becomes very small when U c ≥1.15.In contrast to Fig. 5, the absolute values of slope at the fixed points in Fig. 6 increase with U c .
Hence, the fixed points cannot be an attractor for U c ≥1.In the third panels, Fourier spectra exhibit that except for U c =1, there is only one peak and the predominant frequency increases or the predominant period decreases with increasing U c .This is consistent with Eq. ( 5).Results show that nonlinear behavior disappears when U c >1.
In addition, the amplitude of a peak decreases with increasing U c when U c >1.
Obviously, U c =1 is the critical value of the friction law as displayed in Fig. 4. Fig. 7 exhibits the time variation in V/V max , the phase portrait of V/V max versus U/U max , and Fourier spectrum for four values of η: (a) for η=0.20;(b) for η=0.50;(c) for η=0.87; and (d) for η=0.90 when U c =0.20.In the first panels, the time variation in V/V max exhibits cyclic behavior and the amplitude of V/V max decreases with increasing η.The predominant period of signal only slightly increases with increasing η, because η changes in a small range.In the second panels, the size of phase portrait decreases with increasing U c and there are two fixed points: the first one at V=0 and U=0 and the second one at larger V and larger V. Since the slope values of fixed points are clearly all higher than 1, they are not an attractor.In the third panels, the Fourier spectra exhibit that in addition to the peak related to the predominant frequency, there are numerous peaks associated with higher Ω.This shows nonlinear behavior, mainly caused by nonlinear friction, of the model.The highest peak for case (a) appears at the second frequency.When η<0.9, the amplitude of a peak decreases with increasing η.The frequencies related to the peaks do not change remarkably, because η varies in a small range.Except for case (a), the amplitude of a peak decreases with increasing Ω.The third peak amplitude disappears when η>0.5.The amplitude of a peak becomes very small when Ω>0.25.Except for U c =0.1, the frequencies related to the peaks in Fig. 7 are different from and slightly smaller than those in Fig. 5.Note that when U c <0.55 the simulation results in Fig. 5 are similar to those in Fig. 6.Fig. 8 shows the time variation in V/V max , the phase portrait of V/V max versus U/U max , and Fourier spectrum for four values of η: (a) for η=0.46;(b) for η=0.47;(c) for η=0.98; and (d) for η=0.99 when U c =0.55.When η≤0.47, the time variation in V/V max exhibits cyclic oscillations specified with different main angular frequencies.
When η>0.47 (for example η=0.98 in the figure), in addition to cyclic behavior there is small intermittent slip with shorter periods.This phenomenon also exists when η l <η<η u =0.98.There are unstable (or chaotic) slip when η>η u .Hence, the phase portraits in the second panels display unstable slip at small V and U when η l <η≤η u =0.98.When η=0.99, the velocity becomes an abnormally large negative value at a certain time and the phase portrait also displays unstable or chaotic slip at small V and U.This exhibits unstable slip of the system.In other words, the problem becomes ill-posed in this parameter regime.Since the slope values of fixed points at large V and U are clearly higher than 1, they are not an attractor.Due to the appearance of infinity velocity when η=0.99, the Fourier spectrum is not calculated for η=0.99.The Fourier spectra exhibit that when η<0.47, in addition to the peak related to the predominant frequency, there are numerous peaks associated with higher Ω.This shows nonlinear behavior of the model caused by nonlinear friction.The amplitude of a peak decreases with increasing U c and the peak amplitude decreases with increasing Ω.When η=0.98, the amplitude of the highest peak is much larger than others.For the first three cases, the amplitude of a peak becomes very small when Ω>0.25.The frequencies related to the peaks in Fig. 8 are different from and slightly smaller than those in Fig. 7.
Figs. 9-12 show a variation from stable slip to intermittent slip and then to unstable or chaotic slip when η increases from a smaller value to a larger one for U c =0.6, 0.7, 0.8, and 0.9.The values of η u for U c =0.20-0.95 with a unit difference of 0.05 are given in Table 1.Like Fig. 8, when η≤η l , the time variation in V/V max exhibits only cyclic oscillations specified with different frequencies.When η l <η≤η u , there are small intermittent displacements appear in the cyclic oscillations.Hence, the phase portraits display that unstable slip at small V and U when η l <η≤η u .When η>η u , the velocity becomes an abnormally large negative value at a certain time and the phase portrait also displays unstable or chaotic slip at small V and U.This exhibits unstable slip of the system.In other words, the problem becomes ill-posed in this parameter regime.Due to the appearance of abnormally large negative velocity, the Fourier spectrum is not calculated for η>η u .When η<η l , in addition to the peak related to the predominant frequency, there are numerous peaks related to higher Ω.
This shows nonlinear behavior, mainly caused by nonlinear friction, of the model.The amplitude of a peak decreases with increasing U c and the amplitude of a peak decreases with increasing Ω.For the first three cases, the amplitude of a peak becomes very small when Ω>0.25.Figs.7-12 show that the frequencies related to the peaks slightly decrease with increasing U c and the decreasing rate decreases with increasing U c .In other word, the frequencies related to the peaks for large U c are almost similar.The number of higher peaks and the amplitudes of peaks at higher Ω both decrease with increasing η.This indicates that viscosity makes a stronger effect on higher-frequency waves than on lower ones, and the effect increases with η.
Fig. 13 exhibits the data points of η l (with a solid square) and that of η u (with a solid circle) for several values U c .The values of η l and η u for several values of U c are given in Table 1.The figure exhibits a stable regime when η≤η l , an intermittency (or transition) regime when η l <η≤η u , and unstable (or chaotic) regime when η>η u .
When U c <0.55, there is no η l , in other word, unstable slip does not exist.Clearly, η l , η u , and their difference η u -η l all decrease with increasing U c .This means that it is easier to yield unstable slip for larger U c than for smaller U c .Since smaller U c is associated to larger γ of decreasing rate of friction force with slip, it is easier to yield unstable slip from smaller γ than from larger γ.
Huang and Turcotte (1990Turcotte ( , 1992) ) observed intermittent phases in the displacements based on a two-body model.In other word, some major events are proceeded by numerous small events.Those small events could be foreshocks.They also claimed that earthquakes are an example of deterministic chaos.Ryabov and Ito (2001) also found intermittent phase transitions in a two-dimensional one-body model with velocity-weakening friction.Their simulations exhibit that intermittent phases appear before large ruptures.From numerical simulations of earthquake ruptures using a one-body model with a rate-and state-friction law, Erickson et al. (2008) found that the system undergoes a Hopf bifurcation to a periodic orbit.This periodic orbit then undergoes a period doubling cascade into a strange attractor, recognized as broadband noise in the power spectrum.From numerical simulations of earthquake ruptures using a two-body model with a rate-and state-friction law, Abe and Kato (2013) observed various slip patterns, including the periodic recurrence of seismic and aseismic slip events, and several types of chaotic behavior.The system exhibits typical period-doubling sequences for some parameter ranges, and attains chaotic motion.Their results also suggest that the simulated slip behavior is deterministic chaos and time variations of cumulative slip in chaotic slip patterns can be well approximated by a time-predictable model.In some cases, both seismic and aseismic slip events occur at a slider, and aseismic slip events complicate the earthquake recurrence patterns.The present results seem to be comparable with those made by the previous authors, even though viscosity was not included in their studies.This suggests that nonlinear friction and viscosity play the first and second roles, respectively, on the intermittent phases.The intermittent phases could be considered as foreshocks of the mainshock which is associated with the main rupture.Simulation results exhibit that foreshocks happen for some mainshocks and not for others.

Conclusions
In this work, the multistable slip of earthquakes caused by slip-weakening friction and viscosity has been studied based on the normalized equation of motion of a one-degree-of-freedom spring-slider model in the presence of the two factors.The friction is caused by thermal pressurization and decays exponentially with displacement.The major model parameters are the normalized characteristic distance, U c , for friction and the normalized viscosity coefficient, η, between the slider and the background moving plate, which exerts a driving force on the former.Analytic results at small U suggest that there is a solution regime for η and γ (=1/U c ) to make the slider slip steadily.Numerical simulations lead to the time variation in V/V max , the phase portrait of V/V max versus U/U max , and Fourier spectrum.Results show that the time variation in V/V max , obviously depends on U c and η.The effect on the amplitude is stronger from η than from U c .When U c >1, the time variation of V/V max exhibits cyclic oscillations with a single period and the amplitude of V/V max remarkably decreases with increasing U c .When U c <1, slip changes from stable motion to intermittency and then to unstable motion when η increases.For a certain U c , the three regimes are controlled by a lower bound, η l , and an upper bound, η u , of η.
When U c <0.55, η u is absent and thus unstable or chaotic slip does not exist.When U c ≥0.55, the plots of η l and η u versus U c exhibit a stable regime when η≤η l , an intermittency (or transition) regime when η l <η≤η u , and unstable (or chaotic) regime when η>η u .The values of η l , η u , and η u -η l all decrease with increasing U c , thus suggesting that it is easier to yield unstable slip for larger U c than for smaller U c or larger η than for smaller η.The phase portraits of V/V max versus U/U max exhibit that there are two fixed points: The first one at large V/V max and large U/U max is not an attractor for all cases in study; while the second one at small V/V max and small U/U max can be an attractor for some values of U c and η.When U c <1, the Fourier spectra calculated from simulation velocity waveforms exhibit several peaks rather than one, thus suggesting the existence of nonlinear behavior of the system.When U c >1, the related Fourier spectra show only one peak, thus suggesting linear behavior of the system.
of Nonlinear Processes in Geophysics), Prof. J.G.  1. One-body spring-slider model.In the figure, u, K, η, F D , N, and F denote, respectively, the displacement, the spring constant, the viscosity coefficient, the driving force, the normal force, and the frictional force.

V
/V max versus U/U max , and Fourier spectrum based on different values of model parameters are displayed in Figs.5-12.In the figures, V max and U max are, respectively, the maximum velocity and displacement for case (a) of each figure, because the maximum values of U and V decrease from case (a) to case (d) in this study.First, the cases excluding viscosity, i.e., η=0, are explored.Fig. 5 is numerically made for four values of U c : (a) for U c =0.1; (b) for U c =0.4; (c) for U c =0.7; and (d) for U c =0.9 when η=0. Fig. 6 is numerically made for four values of U c : (a) for U c =1.00; (b) for U c =1.01; (c) for U c =1.15; and (d) for U c =2.00 when η=0.A comparison between Fig. 5 and Fig. 6 suggests that U c =1 is a transition value of the friction law between two modes of slip as displayed in Fig. 4.Only U c <1 is considered below.Secondly, the cases including both friction and viscosity are studied.Fig. 7 is numerically made for four values of η: (a) for η=0.20;(b) for η=0.50;(c) for η=0.87;

Fig. 8
Fig. 8 is numerically made for four values of η: (a) for η=0.46;(b) for η=0.47;(c) for η=0.98; and (d)  for η=0.99 when U c =0.55.The Fourier spectrum is not calculated for case (d), because the velocity becomes an abnormally large negative value at a certain time and the phase portrait also displays unstable or chaotic slip at small V and U.This exhibits unstable slip of the system.In other words, the problem becomes ill-posed in this parameter regime.The time variation in V/V max exhibits cyclic oscillations specified with three main frequencies when η<η l =0.47.There is intermittency slip with shorter periods when η l <η<η u =0.98.There are unstable slip when η>η u .This phenomenon holds also when 0.55<U c <1. Four examples for η varying from η<η u to η>η u for different values of U c are displayed in Figs.9-12.Fig. 9 is made for four values of η: (a) for η=0.39;(b) for η=0.83;(c) for η=0.84; and (d) for η=0.85 when U c =0.6.Fig. 10 is made for four Fig.5exhibits the time variation in V/V max , the phase portrait of V/V max versus

Fig. 6
Fig. 6 exhibits the time variation in V/V max , the phase portrait of V/V max versus U/U max , and Fourier spectrum for four values of U c : (a) for U c =1.00; (b) for U c =1.01;

Figure 2 .
Figure 2. The two types of viscous materials: (a) for the Kelvin-Voigt model and (b) for the Maxwell model.(κ=spring constant and υ=coefficient of viscosity)

Figure 4 .
Figure 4.The plot of η versus 1/U c exhibits the phase portrait and root structure of the system.The solid line displays the function: D(η,1/U c )=η 2 -4(1-1/U c )=0.The solid circle, open circle, and solid square represent, respectively, a stable inflected node with D=0, a stable spiral with D<0, and a stable node with D>0.

Figure 6 .
Figure 6.The time variation in V/V max , the phase portrait of V/V max versus U/U max , and power spectrum for four values of U c : (a) for U c =1.00; (b) for U c =1.01; (c) for U c =1.15; and (d) for U =2.00 for the TP law of F(U)=exp(-U/U c ) when η=0.

Figure 7 .
Figure 7.The time variation in V/V max , the phase portrait of V/V max versus U/U max , and power spectrum for four values of η: (a) for η=0.20;(b) for η=0.50;(c) for η=0.87; and (d) for η=0.90 when U c =0.20 for the TP law of F(U)=exp(-U/U c ).
: υ=υ o exp[(E o +pV a /RT)] where υ o is the largest viscosity at low ambient T of an area, E o is the activation energy per mole, p is the pressure, V a is the activation volume per mole, and R is the universal gas constant 1/2 , τ=ω o t, U=u/D o , U c =u c /D o , andΓ D =F D /K.This gives du/dt=[F o /(mK) 1/2 ] dU/dτ,d 2 u/dt 2 =(F o /mK)d 2 U/dτ 2 .The driving velocity becomes V p =v p /D o ω o Hence, the normalized acceleration and velocity are, respectively, A=d 2 U/dτ 2 and V=dU/dτ.The phase ωt is replaced by Ωτ, where Spray, and one anonymous reviewer for their valuable comments and suggestions to substantially improve this article.The study was financially supported by Academia Sinica, the Ministry of Science and Technology (Grant No.: MOST-105-2116-M-001-007), and the Central Weather Bureau (Grant No.: MOTC-CWB-106-E-02).Table1.Values of η l , η u , and V max for various U c .