Deriving the slip front propagation velocity with the slip- and slip-velocity-dependent friction laws via the use of the linear marginal stability hypothesis

The slip front propagation (SFP) is an important phenomenon in several scientific and technological fields, including frictional sliding on the interface between two media and crack tip propagation in a medium. From a theoretical viewpoint, researchers have investigated SFP using both discrete Lan93 ; Amu1 ; Amu2 ; Tro and continuum Rad ; Suz19 models. Furthermore, laboratory experiments have also been performed to investigate SFP Ben . Geophysical studies have also contributed to the understanding of the SFP behavior because fault tip propagation can be modeled as SFP. The SFP velocity has attracted considerable research interest. For example, ordinary and slow earthquakes are widely known to exist Oba ; Kat , and they have different SFP velocities; the SFP velocities for the former category of the earthquake are much higher than those of the latter. The origin of this difference remains an important open research question.

The friction law plays an important role in the understanding of SFP. For example, a friction stress depends on the state variable in the rate-and-state dependent friction model Amu2 ; Bar13 ; Rad . Bar-Sinai et al. Bar13 investigated the SFP velocity using a semi-infinite continuum model and the rate-and-state friction law. Furthermore, we note that the friction stress can also depend on quantities such as the slip and slip velocity. For example, Myers and Langer Mye employed a spring-block model with a slip-velocity-weakening law to study the SFP velocity. The rate-and-state dependent friction law with the shear stiffness of the contacts can be interpreted as the slip-strengthening behavior Bar12 . Moreover, the friction law depending on the slip and slip velocity is realized for elasto- and visco-dampers Hib03 ; Hib05 . The slip- and slip-velocity- dependences of friction stress is important in geophysical studies Ika ; Ito .

Additionally, the viscosity of the medium is important in understanding the SFP behavior Ots ; Mye ; Suz19 . Our previous study Suz19 analyzed an infinitely long viscoelastic block on a rigid substrate and achieved the systematic understanding of the SFP velocity in terms of the viscoelasticity and the friction law depending on the slip velocity in a quadratic manner. Despite these studies, a thorough investigation of the SFP velocity with viscoelastic medium and a friction law that depends on the slip and slip velocity has not yet been undertaken.

Spontaneous SFP velocity has been widely understood by regarding the phenomenon of SFP as the evolution of a solution of the governing equation from a stable state into an unstable state and employing the linear marginal stability hypothesis (LMSH) Lan91 ; Mye ; Suz19 . The LMSH has been widely used to investigate the dynamics of fronts or domain walls that propagate spontaneously into an unstable state in a model with a nonlinear governing equation Saa1 ; Saa2 . This approach has been used to obtain the solidification front speed Dee , the chemical reaction front speed Saa1 ; Saa2 , and the slip front velocity between blocks and substrates Lan91 ; Suz19 ; Mye . This hypothesis asserts that even if the governing equations are nonlinear, a linearized model yields sufficiently accurate front behavior, including the correct front propagating velocity. This hypothesis has been applied to the slip front behavior via the linearization of the friction law; the friction laws adopted in such studies have been limited, e.g., the slip-velocity-weakening Lan91 ; the systematic treatment of various types of friction law has not been performed.

This paper is organized as follows. A brief introduction to the LMSH is provided in Sec. II.1, and the model setup is defined in Sec. II.2. Analytical treatment of the problem is presented in Sec. III; in this section several SFP velocities are obtained based on the LMSH, and a discussion of their physical implications are presented. The numerical treatment and implications for slip behavior are discussed in Sec. IV, and the treatment of the LMSH is justified there. A summary and discussion are presented in Sec. V.

We now linearize the terms $F(u)-\tau_{\mathrm{v}}(\dot{u})-\tau_{\mathrm{s}}(u)$ from Eq. (6). We take

where $C_{1}$ and $C_{2}$ are constants. To undertake the linearization (7), we assume that the system is in a critical state before the end-loading in the EUL model, which leads to

Analyses without the critical state assumption have been performed in previous studies Mye ; Bar12 . However, those treatments only provide approximations for the SFP velocities. With the assumption of the critical state, an analytical treatment is possible here.

We note that if the slip direction reverses ($u$ turns to be negative) with increasing time, Eq. (7) cannot be applied. The reason for this is twofold: firstly, the friction stress, $\tau_{\mathrm{v}}(\dot{u})+\tau_{\mathrm{s}}(u)$, changes its sign if slip reversal occurs, and thus, the linearization of Eq. (7) cannot be employed. Secondly, if the slip reversal occurs, the absolute value of the stress acting on the slip-reversal point will exceed the maximum static friction stress. This slip criterion is not considered in the linearization of Eq. (7). We do not consider the slip reversal in the following analytical treatment, as is done in seismology Suz06 .

Using the linearization of Eq. (7), we rewrite the governing equation [Eq. (6)],

We note positive and negative values for $C_{1}$ and $C_{2}$ are permitted.

We assume the plane wave solution for Eq. (9), i.e., $u\sim\exp(\pm i(kx-\omega t))$. The upper (lower) sign describes the intruding (extruding) front propagation, as noted in Sec. II.1. Based on this assumption and using Eq. (9), the dispersion relation is found to be

Considering the real and imaginary parts of the dispersion relation in Eq. (10), we obtain,

respectively. We first note that the solutions $k_{r}\neq 0$ and $\omega_{r}\neq 0$ exist. In the case of $C_{1}$ and $C_{2}$ being positive and $k_{i}=\omega_{i}=0$, the solution set $(k_{r},\omega_{r},k_{i},\omega_{i})=(\sqrt{C_{1}},\sqrt{C_{1}+C_{2}},0,0)$ exists. The solutions $k_{r}\neq 0$ and $\omega_{r}\neq 0$ describe an oscillating solution. The slip reversal occurs with these solutions; such solutions cannot be treated in a framework based on Eq. (9). Thus, this solution set is excluded from the following discussion. We therefore assume $k_{r}=\omega_{r}=0$ in the analytical treatment. With such an assumption, considering Eq. (11), we have

Combining this equation, the expression for growth stability [Eq. (2)] and condition for propagation stability [Eq. (3)], we obtain a cubic equation for $\omega_{i}$,

and a relationship between $k_{i}$ and $\omega_{i}$:

We begin by categorizing the solutions of Eq. (14) in terms of the numbers of positive and negative real solutions. Combining the results with Eq. (15), we are able to investigate the numbers of SFP velocities that the above equations yield.

We start by defining $f_{\mathrm{in}}(\omega_{i})$ as being equal to the left-hand side of Eq. (14) with the upper sign; such a definition permits us to study the intruding front velocity. It is noted that the discussion performed in this section is also directly applicable to the extruding front propagation. To confirm this statement, we define $f_{\mathrm{ex}}(\omega_{i})$ as the left-hand side of Eq. (14) with the lower sign. We then have the relationship

Thus, we see that the number of the real solutions for $f_{\mathrm{ex}}(\omega_{i})=0$ is exactly the same as the number of real solutions for $f_{\mathrm{in}}(\omega_{i})=0$, and the positive (negative) solutions for the former correspond to the negative (positive) solutions for the latter. We therefore investigate only the numbers of real solutions for $f_{\mathrm{in}}(\omega_{i})$ in this work.

To categorize the numbers of the real solutions for $f_{\mathrm{in}}(\omega_{i})=0$ in terms of $C_{1}$ and $C_{2}$, we require an expression for $D$, $f_{\mathrm{in}}(0)$, $f_{\mathrm{in}}(\omega_{\pm}^{\mathrm{in}})$ and $\omega_{\pm}^{\mathrm{in}}$, where

and $\omega_{\pm}^{\mathrm{in}}$ is defined as the solution of

which leads to the solution

(double signs in the same order). We observe that if $D$ is positive $\omega^{\mathrm{in}}_{\pm}$ is real. Additionally, we can obtain

and

We define $\omega_{1}^{\mathrm{in}}$, $\omega_{2}^{\mathrm{in}}$, and $\omega_{3}^{\mathrm{in}}$ as the solutions for $f_{\mathrm{in}}(\omega_{i})=0$; $\omega_{1}^{\mathrm{in}}$ is defined to always take real values, as described in Appendix A. We number the various cases based on the diagram shown in Fig. 2(a). In cases 1, 2, 3, 4, 7, and 8, only $\omega_{1}^{\mathrm{in}}$ is real, and $\omega_{1}^{\mathrm{in}}$ is positive in cases 1, 3, and 4, whereas it is negative in cases 2, 7, and 8. Though the numbers of the positive (negative) solutions are the same for the cases 1, 3, and 4 (2, 7, and 8), we differentiate the cases because the numbers of the real positive solutions for $k_{i}$, which is one or zero, may be different in each case. The solutions $\omega_{1}^{\mathrm{in}}$, $\omega_{2}^{\mathrm{in}}$, and $\omega_{3}^{\mathrm{in}}$ are real in cases 5, 6, 9, and 10. All the solutions are positive in case 5, two are positive and one is negative in case 9, two are negative and one is positive in case 6, and all the solutions are negative in case 10.

We now obtain the analytical forms of the boundaries dividing cases 1–10 in $C_{1}-C_{2}$ phase space. These equations are obtained based on the diagram Fig. 2(a) as follows:

First, using Eqs. (17) and (22), we obtain a boundary referred to as boundary A, which reduces to the parabolic form

From Fig. 2, boundary A represents a boundary between a region in which solutions are characterized as cases 1–2 and a region containing solutions that fall into cases 3–10 in the $C_{1}-C_{2}$ space. Next, from Eqs. (20) and (23), we can conclude that $C_{2}=0$ represents a boundary, which we will call boundary B. Boundary B divides regions of solutions that are represented as case 1 from case 2, and the solutions in cases 3–6 from cases 7–10. Further, we consider the cubic equation for $C_{2}$ from Eqs. (21), (24), and (25):

We can factorize this equation to obtain,

which we can evaluate to find that

are boundaries. The former will be referred to as boundary C, and the latter will be referred to as boundary $\mathrm{D}_{\pm}$, respectively (double signs in the same order). Boundary C is a boundary for regions where solutions are of case 3, case 4, and cases 5–6, and $\mathrm{D}_{\pm}$ divides the regions in which solutions in case 7, case 8, and cases 9–10 exist. Additionally, we can confirm analytically that the boundaries A, C, and $\mathrm{D}_{+}$ all pass through the point $(C_{1},C_{2})=(-1.5,0.5)$, and that the boundary C is tangential to the boundaries A and $\mathrm{D}_{+}$ at this point.

Using the boundaries $\mathrm{A}-\mathrm{D}_{\pm}$, we obtain the phase space that categorizes the behavior of the solutions for $f_{\mathrm{in}}(\omega_{i})=0$ (see Fig. 3). Figure 3 shows that cases 6, 7, 8, and 10 are further divided into two disconnected regions. Actualy, from Eqs. (19), (26), and (27), we can obtain the relationship $C_{2}=C_{1}/3$ for $C_{1}>0$ [Eq. (27)] and for $C_{1}<0$ [Eq. (26)]. This forms a boundary between regions in which solutions for $f_{\mathrm{in}}(\omega_{i})=0$ fall into cases 5 and 6, and cases 9 and 10. However, we can see that the case 5 does not appear in the phase space (see Fig. 3), and case 9 is not adjacent to case 10. Therefore, this straight line does not appear as a boundary in Fig. 3.

The solutions for $f_{\mathrm{in}}(\omega_{i})=0$ have been categorized in the $C_{1}-C_{2}$ space. We now categorize the solutions for $k_{i}$. Notably, the solution for $k_{i}$ is obtained from Eq. (15); thus we see that the solution does not exist when the right-hand side of Eq. (15) is negative, even if real and positive solutions for $\omega_{i}$ exist. Therefore, we see that the curve $k_{i}=0$ is another boundary in the $C_{1}-C_{2}$ space. We refer to this boundary as boundary E, and the analytical form of the boundary E is obtained here. If we substitute $k_{i}=0$ to Eq. (13) with the upper sign, we obtain

where the sign of the second term on the right-hand side is arbitrary. Substituting $k_{i}=0$ in Eq. (15) with the upper sign yields

From the two previous expressions, an equation describing boundary E can be obtained:

It can also be confirmed that boundary C is a tangent to boundary E at the point $(C_{1},C_{2})=(-2,1)$. The phase space is shown in Fig. 4. This figure shows 17 subdivided regions of the cases, and they are named as shown in the figure (see Table 1 for details).

We show the values of $\omega_{j}^{\mathrm{in}}$, $\omega_{j}^{\mathrm{ex}}$, $k_{j}^{\mathrm{in}}$, and $k_{j}^{\mathrm{ex}}$ ($j=1,2,3$) in Fig. 5, where $\omega_{j}^{\mathrm{ex}}$, $k_{j}^{\mathrm{in}}$, and $k_{j}^{\mathrm{ex}}$ are defined in Appendix A. Note that $\omega_{2}^{\mathrm{in}}$ and $k_{2}^{\mathrm{in}}$ do not exist. The values of $k_{3}^{\mathrm{in}}$, $k_{1}^{\mathrm{ex}}$, and $k_{3}^{\mathrm{ex}}$ are zero on boundary E. Note that since Eq. (34) is a necessary and not a sufficient condition for $k_{j}^{\mathrm{in}}=0$ or $k_{j}^{\mathrm{ex}}=0$, non-zero values ($k_{1}^{\mathrm{in}}$ and $k_{2}^{\mathrm{ex}}$) are allowed on bounary E.

Using the values of $\omega_{j}^{\mathrm{in}}$, $\omega_{j}^{\mathrm{ex}}$, $k_{j}^{\mathrm{in}}$, and $k_{j}^{\mathrm{in}}$ corresponding to solutions of $f_{\mathrm{in}}(\omega_{i})=0$, we also investigate the intruding front velocity, $v^{\mathrm{in}}$, and the extruding front velocity, $v^{\mathrm{ex}}$. We define the $j$th intruding and extruding front velocity as $v_{j}^{\mathrm{in}}$ and $v_{j}^{\mathrm{in}}$, respectively, as in Appendix A, and the velocities $v^{\mathrm{in}}_{1}$, $v^{\mathrm{in}}_{3}$ $v^{\mathrm{ex}}_{1}$, $v^{\mathrm{ex}}_{2}$, and $v^{\mathrm{ex}}_{3}$ exist. The regions where $v_{1}^{\mathrm{in}}$, $v_{3}^{\mathrm{in}}$, $v_{1}^{\mathrm{ex}}$, $v_{2}^{\mathrm{ex}}$, and $v_{3}^{\mathrm{ex}}$ exist are shown in Fig. 6(a–e), respectively, in $C_{1}-C_{2}$ space. We note that the region where $v_{3}^{\mathrm{ex}}$ exists is separated into two disconnected regions. We refer to $v_{3}^{\mathrm{ex}}$ which corresponds to values of negative $C_{1}$ as $v_{31}^{\mathrm{ex}}$, and that which corresponds to positive values of $C_{1}$ as $v_{32}^{\mathrm{ex}}$. We have three velocities in cases 6a and 9b ($v_{1}^{\mathrm{in}}$, $v_{2}^{\mathrm{ex}}$, $v_{32}^{\mathrm{ex}}$ for case 6a and $v_{1}^{\mathrm{in}}$, $v_{3}^{\mathrm{in}}$, $v_{1}^{\mathrm{ex}}$ for case 9b). There exist two velocities in cases 9a and 10d ($v_{1}^{\mathrm{in}}$ and $v_{2}^{\mathrm{ex}}$ for 9a and $v_{1}^{\mathrm{ex}}$ and $v_{31}^{\mathrm{ex}}$ for 10d). For cases 1, 4, and 6b, we obtain the single velocity $v_{1}^{\mathrm{in}}$, whereas in cases 2a, 7a, 8b, and 10b, only a single velocity $v_{1}^{\mathrm{ex}}$ is found to exist. Finally, we note that no solutions for $v^{\mathrm{in}}$ or $v^{\mathrm{ex}}$ exist in cases 2b, 7b, 8a, 8c, 10a, and 10c.

The values of the front velocities $v_{1}^{\mathrm{in}}$, $v_{3}^{\mathrm{in}}$, $v_{1}^{\mathrm{ex}}$, $v_{2}^{\mathrm{ex}}$, $v_{31}^{\mathrm{ex}}$, and $v_{32}^{\mathrm{ex}}$ are shown in Fig. 7. It is found that $v_{3}^{\mathrm{in}}$, $v_{1}^{\mathrm{ex}}$, and $v_{31}^{\mathrm{ex}}$ diverge close to boundary E. This divergence is because the values of $k_{3}^{\mathrm{in}}$, $k_{1}^{\mathrm{ex}}$, and $k_{3}^{\mathrm{ex}}$ are equal to zero, while $\omega_{3}^{\mathrm{in}}$, $\omega_{1}^{\mathrm{ex}}$, and $\omega_{3}^{\mathrm{ex}}$ take non-zero values, on the boundary. We can therefore conclude that these velocities describe optical modes and represent unphysical propagations. In addition, we confirmed that the velocities $v_{1}^{\mathrm{in}}$ and $v_{2}^{\mathrm{ex}}$ take the values $\sqrt{1+C_{1}}+\sqrt{C_{1}}$ and $\sqrt{1+C_{1}}-\sqrt{C_{1}}$, respectively, on the $C_{1}$ axis for $C_{1}>0$. These values are consistent with those found in Ref. Suz19 because the slip-dependence of the friction stress was negelected in that work, i.e., $C_{2}=0$. The value of $v_{32}^{\mathrm{ex}}$ is not consistent with $\sqrt{1+C_{1}}-\sqrt{C_{1}}$. We therefore conclude that $v_{1}^{\mathrm{in}}$ and $v_{2}^{\mathrm{ex}}$ represent the physical SFP velocities for the intruding and extruding fronts, respectively. We can also suggest that $v_{1}^{\mathrm{in}}$ and $v_{2}^{\mathrm{ex}}$ are the extension of the intruding and extrding front velocities obtained in Ref. Suz19 to the model with the friction law depending on the slip and slip velocity.

In this section we show the spatiotemporal evolutions of the slip velocity found via numerical computations using the values shown in Table 1. We confirm that the values of $v_{1}^{\mathrm{in}}$ and $v_{2}^{\mathrm{ex}}$ obtained using the LMSH are accurate.

In this work we have considered an infinite 1D viscoelastic block on a substrate subject to end-loading stress and the end- and upper-loading stresses together with a friction law that depends on the slip and slip velocity. The two intruding front velocities and three extruding front velocities were analytically obtained based on the LMSH, which requires only the linearized governing equation. Of these five velocities, three represent the optical modes, and the unique intruding and extruding front velocities were found to exist for physically meaningful SFPs. Via numerical calculations, the intruding front velocity was obtained for the EUL model, and the extruding front velocity was realized for the EL model; these results indicated that the detail of the dependence of the friction stress on the slip and the slip velocity does not affect the front velocity, and that the linearized governing equation is sufficient to obtain accurate results.

Our framework is found to be quantitatively consistent with the results of previous studies. For example, the dependences of $v_{1}^{\mathrm{in}}$ and $v_{2}^{\mathrm{ex}}$ on $C_{1}$ for $C_{2}=0$ are consistent with those of $v^{\mathrm{in}}$ and $v^{\mathrm{ex}}$ in Ref. Suz19 , respectively, as mentioned in Sec. III. Additionally, the form of $v_{1}^{\mathrm{in}}(C_{1})$ for $C_{2}=1$ is consistent with that obtained in Myers and Langer Mye . Though they obtained only an approximate solution, the exact analytical solution has been obtained here. Furthermore, note that their model corresponds to the case $C_{2}=1$, and they changed parameters $\eta$ and $\alpha$, where $\eta$ is the viscosity, and $\alpha$ corresponds to $C_{1}$. Therefore, both their and our models have two parameters. However, our model is superior because we can treat both positive and negative $C_{2}$ in a single framework. Moreover, we were successful in explaining the variation in SFP velocities using only friction laws with constant viscosity.

We also have shown the regions in which the SFP velocities exist in the $C_{1}-C_{2}$ phase space; we also have presented the analytical forms of the boundaries of these regions in the $C_{1}-C_{2}$ phase space. This phase space may be useful to predict, e.g., whether ordinary or slow earthquakes are likely in natural faults, if the friction law is determined.

Since LMSH linearizes the governing equations, it cannot give an estimation of the macroscopic slip profile. Note that the pulse-like slip was observed in Fig. 9, and the step-function-like slip emerged in Fig. 10. These are consistent with previous laboratory experiments Nie or seismological observations Das ; Aag . Nonetheless, predicting these slip behaviors is difficult with the present framework.

Finally, the critical state was assumed for the EUL model. Based on this assumption, we can proceed with the analytical treatment. In reality, the SFP in a non-critical state has been analyzed Mye ; Bar12 , and the SFP velocities have been obtained. They are, however, approximations. Though the current treatment cannot be directly applied to the non-critical case, we anticipate that it will be an important step in analytically treating SFP velocities with the non-critical state.