How does dissipation affect the transition from static to dynamic macroscopic friction?

Description of the transitional process from a static to a dynamic frictional regime is a fundamental problem of modern physics. Previously we developed a model based on the well-known Frenkel-Kontorova model to describe dry macroscopic friction. Here this model has been modified to include the effect of dissipation in derived relations between the kinematic and dynamic parameters of a transition process. The main (somewhat counterintuitive) result is a demonstration that the rupture (i.e. detachment front) velocity of the slip pulse which arises during the transition does not depend on friction. The only parameter (besides the elastic and plastic properties of the medium) controlling the rupture velocity is the shear to normal stress ratio. In contrast to the rupture velocity, the slip velocity does depend on friction. The model we have developed describes these processes over a wide range of rupture and slip velocities (up to 7 orders of magnitude) allowing, in particular, the consideration of seismic events ranging from regular earthquakes, with rupture velocities on the order of a few km/s, to slow slip events, with rupture velocities of a few km/day.


nonlinear di
sipative process.It is generally accepted that friction appears due to interactions between surface asperities.The actual contact area between rough frictional surfaces of stiff materials is less (usually much less) than the nominal surface area, is proportional to the averaged normal stress, and depends on the elasticity and plasticity of the materials in contact [1,2].The normal stress at the tip of an asperity (i.e., at the physical contact area) is equal to the penetration hardness of the material [2].Under static or uniform sliding conditions, friction is usually described by the frictional coefficient, i.e. the proportionality coefficient between tangential and normal stress (classical Amontons-Coulomb law).However, as shown in modern laboratory experiments, friction depends on slip, sliding rate, contact time and normal stress history (see extensive reviews by Marone [3]), Baumberger and Caroli [4] and Dieterich [5]).Recent laboratory experiments [6][7][8][9] also confirm that such a description is not sufficient when parameters describing the physical state of the system (sliding rate, stress, etc.) are not uniform in time and/or space.

Over the past 50 years, various approaches for the modeling of non-uniform frictional processes have been de eloped.Two types of models are the most common, i.e. mass-spring models [10][11][12][13][14][15][16][17], and rate-and-state (Dietrich-Ruina) models [18][19][20][21][22][23][24][25][26][27][28].The mass-spring models of the Burridge-Knopoff type describe collective behavior and statistical features of earthquakes and reproduce major empirical laws of observed seismicity, i.e., large earthquake recurrence, the Gutenberg-Richter law, foreshock and aftershock activities, and preseismic quiescence [10,11,13].In general, however, these models are not adequate to describe the dynamics of an individual event.More detailed dynamics, such as the necessary and sufficient conditions for nucleation of individual earthquake events, have instead been formulated in the framework of rate-and-state models [19][20][21].Rate-and-state models have been used successfully to describe regular earthquakes [18][19][20][21][22][23], slow slip events [24,25], and fault dynamics [e.g.22,23,26].They are capable of incorporating such phenomena as frictional dilatancy [29][30][31], compaction of brittle materials [29,32] and microscopic elasticity [33].Although these models are based on laboratory experiments and include some measurable laboratory parameters, such as the characteristic slip distance, they also include unknown parameters which can be adjusted to fit field or laboratory observations.The same is true for the mass-spring models.Ultimately, a physical model with no adjustable parameters is most desirable.

The Frenkel-Kontorova (FK) model [34] provides a promising point of departure for a more predictive model.It has been idely used to describe micro-and nano-scopic friction (e.g.[35,36] and references therein).Recently, we have developed a FK-type model which describes macroscopic friction.The advantages of this model are: 1) it is an intrinsically dynamical model, rooted in the Newtonian equations of motions; 2) parameters used in the model have explicit and unambiguous physical correlates; 3) it describes frictional processes over a wide range of conditions, from very fast processes such as regular earthquakes down to very slow processes such as creep, silent, and slow earthquakes [37][38][39].The observed nonlinear dynamics of frictional processes is incorporated in the standard linear mass-spring models by introducing ad hoc nonlinear relations between various model parameters (e.g., the introduction of a nonlinear spring constant or a nonlinear relation between friction and slip velocity).By contrast, the FK model is inherently nonlinear.

The motivation for using this model to describe dry macroscopic friction derives from the similarity between plasticity and ry friction, both on laboratory and geophysical scales.Sporadic local motions of the Earth's crust along faults, occurring due to earthquakes and various creeping events, are similar to the processes of plastic deformation in crystals resulting from movement of edge dislocations by the localized shift of crystalline planes.Of particular relevance is the fact that the external stress initiating plasticity is only a small fraction of the stress necessary for the uniform relative displacement of planes of crystal atoms.Similarly, laboratory friction experiments have shown that the critical shear force needed to initiate a macroscopic slip pulse between frictional surfaces is usually much less than that predicted by theory [7].

Motivated by these similarities, we have proposed a novel model [37,38] in which sliding occurs in much the same way as in plas icity, i.e., due to movement of a certain type of defect, a "macroscopic dislocation," which requires much less shear stress than uniform displacement of frictional surfaces.A dislocation is a static configuration of accumulated stress between two surfaces due to the elastic shift of asperities on one surface relative to the other.The dislocation is confined to a specific region of the surface but can be displaced along the surface.As we will see, the sliding motion of the two surfaces occurs due to the movement or propagation of dislocations, somewhat analogous to caterpillar motion.A macroscopic slip is, in fact, the result of a multitude of dislocations propagating along the surface.

In the continuum limit, our model is described by the sine-Gordon (SG) equation, one of the fully integrable nonlinear equations of mathematical physics.This equation has been thoroughly investigated due to its exceptional importance and universality [40][41][42][43].The mathematical apparatus which has been developed is fully applicable to the problems considered here.In the framework of our model, all variables, whether or not directly measurable, are connected by transcendental analytical relations, allowing a clear analysis of dependencies, e.g., rupture velocity as a function of accumulated shear stress [37,38].Algebraic relations have been obtained between kinematic parameters (such as slip velocity and rupture velocities) and dynamic parameters (such as shear stress, normal stress, and stress drop).However these formulae (analytical solutions) have been derived neglecting friction.Here, we introduce a dissipative term into the SG equation, which requires a numerical solution of the problem.We show that some of our previous results, such as the relation between rupture and shear stress, remain valid under the influence of dissipative processes, but some of them, such as the relation between slip velocity and shear stress, need to be modified.

In the next section we describe the basics of the model, followed by the Results, Discussion, and Conclusion.


Model

An overview of the odel is provided to establish the context for the results which follow.A more detailed description may be foun

in our
previous articles [37,38].


Model derivation

Asperities on each frictional surface are idealized as uniform sinusoidal surfaces, illustrated in Figure 1.We will cons

er asperities on
ne of the frictional surfaces as forming a linear chain of balls of mass M, each ball interacting with its nearest neighbors via spring forces of stiffness K b .These provide the forces of elastic deformation for the shift of an asperity from its equilibrium position.The asperities on the opposite frictional surface are regarded as forming a rigid substrate which interacts with the masses M via a sinusoidal restoring force.The physical correlate of this force is the horizontal component of the normal force exerted by the lower surface on upper plate asperities displaced as in figure 1a.Application of this model to describe the slip dynamics yields the FK model, where we have also included an explicit frictional force f i on the i th asperity:
  2 11 2 2 ( 2 ) sin ( , , ), , ii b i i i d i i uu M K u xt u u F u F f x t t b t           (1)
where u i is the shift of ball (a he amplitude of the periodic restoring force and F is the external (or driving) force.In this model, only the interaction between nearest-neighbor asperities is considered.This approximation is motivated by the fact that a disturbance of the stress field around an asperity decreases with distance as r -3 , thus the interaction between nearest neighbor asperities is at least 8 times larger than the interaction between the next closest asperities twice as far away.Following the same procedure used to describe plasticity [42,44,45], we express the coefficients of equation ( 1) through the parameters of the material and the frictional surfaces.For a volume density ρ we find: 3 ,
Mb   2 , (1) b b K     2 , 2 d b F   
where μ is the shear modulus and ν the Poisson ratio.

Then equation ( 1) can be written in the conti ) ( ) , ( / ) ( / ) u b u b u A A F f tc b x b b b           (2 ionless parameter A is equal to
2 / 1 )] 2 / ) 1 [(  
. Note that in deriving the FK model to describe plast s essentially the ratio of the amplitude of two forces: one is the force amplitude between an atom and the substrate layer and the other is the force amplitude between neighboring atoms at the top layer.To describe the respective coefficients for the situation where slip occurs between two external surfaces in contact, we will use equation ( 2) with one significant change: we shall treat the parameter A phenomenologically, using the result for a crystal as a guide.So we assume that A likewise depends on the ratio of two relevant forces.The force amplitude experienced by an asperity due to neighboring asperities along the slip direction is exactly the same as it was for the case of plasticity.But the force amplitude between asperity and substrate is different and depends on the normal stress Σ N .Indeed, when Σ N = 0 the force is zero, since there is no interaction between asperities and a substrate.On the other hand, when Σ N reaches the penetration hardness p  , the interface between the two blocks disappears and the corresponding force amplitude is essentially the same as in the case of plasticity.So we regard A as a function of the ratio of Σ N to :
p  ( / ). Np Af   The simplest choice is 1/2 ((1 ) / 2) NN pp A      
.
ces interpenetrate and it is the ratio between actual and nominal contact areas [1,2].Equation ( 2) in dimensionless form is


Uniform sliding motion

Equation (3) in the absence of external and frictional forces is the w e of a classical wave solution traveling to the right with wave velocity U (in units of nits of A/b), in which  is a function of  −  .Define  =  −   and  =  .The basic solutions are phonons, breathers, and kinks (a particular class of solitons) [41,42].Only the soliton solutions represent the wave propagation necessary to model frictional dynamics.These solutions are characterized by |U|<1.Integrating equation (3) fo its derivatives in terms of the elliptic Jacobi functions, cn and dn [41]:
𝑢 = arcsin[±𝑐𝑛(𝛽 𝜉)], 𝜎 𝑠 = 2𝛽 𝑑𝑛(𝛽 𝜉), 𝑤 = 𝑈𝜎 𝑠 ,(4)𝑘 = 2𝜋𝑁 = 𝜋𝛽 𝐾(𝑚) ,
where () is the complete elliptic integral of the first kind of modulus 
and N is the density of kinks in units of A/b.Solution (4) describes an infinite sequence of interacting kinks (solitons) of one sign, which are periodic in space and time.In the context of our model, a soliton is a dislocation.It is also useful to introduce three dimensionless variables averaged over an oscillation period
, , . 2 2 2 2 s S d wd UN d W k E                 (5)
These variables correspond to the measurable parameters of slip velocity, stress and strain.The parameters of a dislocation (amplitude of stress 0 s  and strain 0  associated with the presence of the dislocation) are:
   / 0 A s  , ) 2 /( 0   A  .
Let us describe a scenario of frictional processes in terms of solutions ( 4) and ( 5).The macroscopic dislocations are nucleated on the surface by an applied shear stress in the presence of asperities.As in crystals, the mobility of a macroscopic dislocation over the frictional surface is much larger than the mobility of the whole surface, since the displacement of a dislocation (a pre-stressed area) requires less external stress.So the relative sliding of two bodies occur ment of dislocations.The passage of a dislocation through a particular point on the sliding surface shifts the contacted bodies locally by a typical distance b.Such a dislocation may propagate with any velocity U ranging from 0 to c, and the a

rage veloc
ty of sliding, i.e. the observable slip rate W, is proportional to the dislocation velocity U and dislocation density N. The parameters of a dislocation (stress amplitude and pulse width) are entirely defined by the material parameters and the normal stress and do not depend on process parameters such as dislocation density and slip rate.The characteristic width of a dislocation is
A b D / 2 
. Usually, the dislocation width is much larger than the typical distance between asperities ranging from 10 2 to 10 4 of the asperity size.Figure 2 schematically illustrates the frictional processes via movement of macroscopic dislocations. .W is the slip velocity w averaged over an oscillation period in time.By contrast, the w velocity field is spatially non-uninform at smaller scales, i.e., on a scale comparable with a dislocation width (panel (b)).The spatial and temporal averages of w in (b) give the uniform field in (a).The dislocation size is

ually much larger than the
ypical size of an asperity.The relative movement on the frictional surfaces at even smaller scales (asperity size) could be larger than average if it is at the center or peak amplitude of the dislocation (panel (c)) or smaller than average if it is in between two dislocations (panel (d)).The values in the panels reflect experiments described in [6][7][8][9].


Non-uniform sliding motion

Solutions ( 4) and ( 5) can be used to describe the frictional process in the macroscopically uniform case.However, in the non-uniform case the macroscopic parameters such as slip velocity and stress are functions of time and space.In typical cases, the frictional surface includes many dislocations.A numerical solution for the exact positions of a multitude of microscopic dislocations is too computationally demanding to be practicable.Moreover, they cannot be measured and are therefore of little value except as they contribute on the average to the measured macroscopic dislocation.We therefore apply a technique which simplifies equation (3) at the expense of losing non-essential details of the process such as the exact position of each dislocation.Witham [46] developed a procedure for constructing a system of equations describing the dynamics of averaged variables (i.e., the variables in equation 5) that derive from the original dynamical equations (i.e., equation ( 3)).For strict applicability, the average values should vary slowly in time and space.
r the present application, the variables of interest can all be expressed in terms of the two independent variables U and m.Applying this procedure to the homogeneous (  −  = 0) SG equation yields the system of coupled equations [45]:
22 22 11 1 0, 1 2 1 2 1 0, 1 2 1 2 t t x x t t x x UU U m U m U m U m UU U m U m U mm U mm             (6) where , K E   m m   1 1
, and is the complete elliptic integral o nces [43,45].However, system (6) does not include friction.

Whitham [46] also described a formal procedure for including

dissipat
ve term.Applying this procedure we obtain the Whitham equations for equation ( 3):
E m ( ) 2 1/2 2 2 22 11 41 [ ] ( , ), [m(1 )] 1 2 1 2 1 0, 1 2 1 2 t t x x t t x x K U U U m U m D m U U U m U m UU U m U m U mm U mm                 (7) where 1 ( , ) [ (u, u )] 2 t D m U F f d   
. This system of equations does not have analytical solutions (in contrast to system (6)) and must be solved numerically.


Results

We first we describe the results derived from the analytical solutions of equations ( 6).We then

ompute solutions of equations (7), wh
ch includes a dissipative term, and compare the results in order to characterize the effects of friction.Finally, we consider how the initial stress distribution affects the solution.Mathematica was used to obtain numerical solutions of (7).The initial val ed through the definition of Σ  from equations ( 4) and (5).


Analytical solution (no dissipation)

To model the transition from the static to the dynamic regime of the frictional process, we consider the following idealized problem.Suppose the point x = 0 divides the areas of stressed (x < 0) and unstressed (x > 0) material, modeled as a step-function
( 0, 0) , s t x       0 ) 0 , 0 (        x t s , ( 0, ) 0 U t x 
 (solid line in Figure 3).With these initial conditions, we have
Σ ± = 𝜋 √ 𝑚 ± 𝐾(𝑚 ± )
. This setting models a stress accumulation in front of an obstacle placed at position x = 0 (e.g., a large asperity).The effect of alternative initial shear stress profiles is a topic for further investigation.Assume that at time t = 0 the external shear stress reaches the value necessary to overcome this obstacle.Then the dynamics of the transition (solution of system (6)) is described by the following formulae [38,45]:
( ) ( ), , , s x V m U k W kU t Km                     (8) 1 , , ,1G G V V V G G                         (9)
where
) /( ) ( 1 1 m K E m K E G     , m m / ) 1 ( 2 1    , () G G m  and() m    
, and m - is defined by the transcendental relation
() K m m     
. The variable V is the nonlinear group velocity of the wave solution in units of c.Along a line x/t = V = constant in the x-t plane, all variables are constant.The solution is represented by a region expanding in time and bounded between the lines x/t = V(m= 0) = V − = −1 and x/t = V(m=1) = V + .Note that inside the expanding region all variables are functions of time and position.The indices + and − designate the leading and trailing edges of the slip pulse, respectively.Thus the solutions (8 and 9) describe the dynamics of a slip pulse with two rupture fronts (or detachment fronts, in terms of equation ( 3)) propagating in opposite directions with different velocities V − and V + .Figure 4 shows the ratio of V − to V + as a function of V − .One can see that for small rupture velocities the ratio is large, i.e. the pulse propagates practically in one direction.Figure 5 shows the result of the calculation of rupture velocity V -as a function of dimensionless initial shear stress   .Note a strong dependence of velocity on stress.Indeed, changes of shear stress by an order of magnitude lead to velocity changes by seven orders.The velocity of dislocation movement U(x,t) (formulae ( 8)) ranges from zero at the pulse trailing edge to the value V + at the pulse leading edge.The movement of dislocations is accompanied by slip with velocity W(x,t) (formul

bottom panel).The reason f
r this discrepancy arises because dissipation was not taken into account (see next subsection).


Solution with dissipation

Let us consider the same problem as above, i.e. the transition from the static to the dynamic regime, described now by equations (7).Note the small difference in initial conditions (dotted line at Figure 3) compared to the case in section 3.1 (solid line at Figure 3
): 1)      but 0  and
2) the stress distribution is a smooth function of x (not a step-function).These changes were necessary in order to obtain a numerically stable solution.For a simple velocity-dependent frictional term,
t f Cu  and supposing that 0, F  we find ( , ) (m, U), D m U CW 
where C is the dissipation coefficient and W the slip velocity as defined by the formulae (8). Figure 7 depicts the spatial and temporal distribution of a slip pulse velocity for the various values of coefficient C. One can see that the amplitude of the slip velocity decreases with an increase in C . Figure 8 shows the slip velocities at time 2π for various values of C. To examine how friction (the dissipative term) affects the rupture velocity, we normalize W to make the maximum slip velocity the same for all cases.Then we see that both velocities V -and V + (left and right sides of slip pulse) are practically unchanged with increasing coefficient C for all cases considered, although the slip velocity amplitude changes by a factor greater than five.Note that this is true for any value of the i

tial stress   provided    
 .The spatial distribution of shear stress at time 2π is depicted in Figure 9.One can see that the larger dissipation, the smother is the transition from the stressed to the less stressed region.The result of calculating the dependence of slip velocity on the dissipative term for various initial shear stress values is shown in Figure 10.One can see that dissipation (dynamic friction) can essentially reduce the slip velocity.This velocity reduction depends weakly on the v lue of the initial shear stress   .


Influence of initial stress Σ +

We have shown in section 3.1 how the rupture velocities and slip velocity depend on initial stress   .How do the

 and ΔΣ=Σ
--Σ + values affect the shape and parameters of a slip pulse? Figure 11 depicts the spatial and temporal distribution of a slip pulse for various values of   .One can see that the shape of the pulse and its velocity do not change significantly unless the value of Σ + ≈ Σ − .. Analysis show that   does not affect the fast rupture velocity ;

V  however the slow rupture velocity      Regular and slow earthquakes in the Earth's crust are typical transitional processes.Shear stress concentrated along a fault is relaxed due to various types of earthquakes and creep.Observations of far-and near-source ground motion allow reconstruction of kinematic parameters of regular earthquakes, i.e. rupture and slip velocities [47][48][49][50][51]. Dynamic parameters could be found by the modelling of processes using additional information and/or assumptions [51].According to our model, the rupture velocity V -is explicitly defined by the initial stress Σ -and the elastic parameters of the medium, and doesn't depend on the dissipation term and the initial stress Σ + or stress difference ΔΣ.Let us consider, as an example, the 2004 M=6 Parkfield earthquake.According to [51], the rupture velocity was about dim 3.0 V   km/s and the slip velocity at the hypocenter area was about dim 0.5 W  m/s.Taking the P-wave velocity in the Parkfield area to be c l = 6 km/s and ν = 0. .Now we can find the dimensionless slip velocity (Figure 6, upper panel, red dashed lines) to be W=0.1.Recall that the velocity in Figure 6 is calculated for C=0.To find the velocity for other values of C we need to use Figure 10.Thus, for the case considered, the dimensionless velocity is W=0.083, 0.05, and 0.02 for C=0.1, 0.5 and 2.0, respectively.Now we can calculate dim / (cW) 0.0029 AW   , 0.0035, 0.0058 and 0.0145 for the cases C=0.0, 0.1, 0.5 and 2.0, respectively.Supposing that µ = 30 GPa we can calculate the initial stress to be dim
V  increases with increasing   or decreasing ΔΣ. When / 1     , V  is practically equal to , V  i./ 29 A       ,
34, 58 and 145 MPa for these respective values of C. The result of a sophisticated dynamic modelling of the Parkfield earthquake [51] yields a stress value of about 31 MPa, which practically coincides with our estimate for the case of negl ation.This leads us to the unexpected conclusion that (at least for some regular earthquakes) the dissipation due to friction may have practically no effect on the slip velocity, W, in addition to having no effect on the rupture velocities V

nd V + .Now
et us consider an example of slow earthquakes, e.g.so-called episodic tremor and slip phenomena [52][53][54].It is known that during these events, a slip pulse slowly propa ates along a subduction fault with an effective rupture velocity dim 10

V   km/day .To estimate the actual value of W we need to know the value of A, which has been estimated by us in a previous work to be 5 A 4 10   [39].Using this value we find the dimensionless velocity to be 7 dim / (cA) 0.5 10
WW     
. This value is one-fifth the value predicted by the model without dissipation.So for the slow event considered, friction has a pronounced influence on slip velocity.The corresponding dissipative coefficient (see Figure 10) is about C=2.0.


Conclusions

Our model for describing the transition from static to dynamic friction has been ext nded here to include dissipation.The significant and novel elements of the model are:

1) The FK model has been modified and adapted to the case of non-lubricant macroscopic friction by introducing a parameter A which is the ratio between the real and nominal area of the frictional surfaces;

2) In the model proposed, sliding occurs due to movement of a certain type of defect, a "macroscopic dislocation" (or area of localized stress), which requires much less shear stress than uniform displacement of frictional surfaces.

3) To describe measureable macroscopic parameters, we use the Whitham modulation equations applied to the SG equation rather than the SG equation itself which applies to details at the scale of individual dislocations.

The source of the transition from static to dynamic friction is the gradient of the shear stress (Figure 3).We show that in transition processes, the rupture velocity of the trailing edge V -, i.e. the velocity of a detachment front propagating through the stressed material, is defined only by the initial stress   and does not depend on friction in the case of a simple linear dissipative term considered here.This counterintuitive conclusion had already been predicted in our previous article [38], which neglects friction in the calculation, by comparing the accumulated elastic energy and dissipative energy released during the transition.Here we verified this prediction by explicitly including friction in the calculation (see Figures 8 and 11).Note that the initial stress (expressed in dimensional units) is proportional to the ratio between shear and normal stress.Thus, the rupture velocity of the trailing edge V -is defined by the ratio between shear and normal stress.This conclusion is consistent with the experimental result that the velocity of the detachment front is defined by this ratio [9].The rupture velocity of the leading edge V + , i.e. the velocity of the detachment front propagating from stressed to le s stressed or unstressed material, is defined by the initial stresses   and   and likewise does not depend on friction.The velocity V -is always larger than V + (Figure 4), i.e., the rupture front always propagates more easily through stressed rather than unstressed or less stressed material.However, in

Figure 1 .
1
Figure 1.Schematic of asperity contact (a) and chain of masses interacting via elastic springs and placed in a periodic potential (substrate) (b).The balls represent asperities.The sine-shaped surface is the opposite plate.




x and t are now in units of b/(2π), b/A and b/(cA), respectively, F and f are the external force and frictional force per unit area in units of / interpreted as the dimensionless strain and the dimensionless slip in units of ) the driving force is the tangential stress, we set it equal to t