Towards a simple, comprehensive model of regular earthquakes and slow slip events, part I: one-dimensional model (original) (raw)

Towards a simple, comprehensive model of regular earthquakes and slow slip events, part I: one-dimensional model

Naum I. Gershenzon*1, Cong Zhou 2{ }^{2}, and Thomas Skinner 1{ }^{1}
1{ }^{1} Physics Department, Wright State University, 3640 Colonel Glenn Highway Dayton, OH 45435
2{ }^{2} The Second Monitoring and Application Center, No. 316 Xiying Road, Xi’an City, Shaanxi
Province, China, 710054
*Corresponding author
E-mail address: naum.gershenzon@wright.edu

Abstract

Seismic events are a prominent manifestation of the earth’s fault dynamics. Since the focal area of these events is inaccessible to direct observation, modeling plays a vital role in the investigation of processes below the earth’s surface. We have developed a model that describes the major characteristics of a rupture, ranging from regular earthquakes (EQs) to slow slip events (SSEs), including episodic tremor and slip (ETS). Previous model predictions, while accurate, are based on a highly idealized initial stress distribution and a simple velocity-dependent expression for friction. The full scope of the model has, therefore, not been fully demonstrated. Further developments, presented here, include more physically realistic treatments of both the initial conditions and friction. Important new model predictions are: (1) The type of a seismic event, i.e. regular EQ or SSE, is determined by the fault strength, the shear to normal stress ratio, and the gradient in the ratio. Quantitative values for these crucial parameters are also obtained here. The critical factor for the appearance of a regular EQ is the abrupt decrease of friction at slip velocities ≥0.1 m/s\geq 0.1 \mathrm{~m} / \mathrm{s} due to flash heating or other mechanisms. If the stress ratio and spatial

gradient of this ratio are large enough to accelerate slip velocity up to ∼0.1 m/s\sim 0.1 \mathrm{~m} / \mathrm{s}, then a regular EQ occurs. For lower values, all types of SSEs are possible; (2) Rupture velocities for regular EQs range from a fraction of the shear wave velocity, csc_{s}, up to the supershear velocity. The range for SSEs is much wider, from 10−7cs10^{-7} c_{s} to 10−1cs10^{-1} c_{s}. The maximum slip velocity for regular EQs is typically on the order of 1 m/s1 \mathrm{~m} / \mathrm{s}. For SSEs, the slip velocity ranges widely from a few cm/\mathrm{cm} / year up to 0.1 m/s0.1 \mathrm{~m} / \mathrm{s}. The magnitude of the stress drop may vary in the range from 1%1 \% to 10%10 \% of the initial shear stress for regular EQs and from 0.1%0.1 \% to 10%10 \% for the SSEs. Generally, the stress drop in SSEs is one to three orders of magnitude less than in EQs; (3) The rupture can expend as a crack-like mode or as a self-healing pulse mode. The type of rupture mode is determined by the stress ratio and its gradient. If stress heterogeneity has a step-like shape, the rupture is always crack-like. If the heterogeneity is localized, the rupture is crack-like for sufficiently large values of the stress and its gradient (as quantified in this work) and is pulse-like for smaller values. The type of heterogeneity also determines the position of maximum slip; (4) After an instability develops, rupture dynamics do not depend on the relative values of the ratestate aa and bb parameters, i.e., it does not matter if frictional sliding is velocity-weakening ( a<ba<b ) or velocity-strengthening ( a>ba>b ). Popular models based on the rate-state frictional law are not consistent with this result.

1. Introduction

The ultimate goal of earthquake research is a predictive theory for seismic events. Towards that end, we have developed a model and investigated its relevance for studying the dynamics of earth’s crustal faults [1-8]. The model is derived solely from physical principles, has no adjustable parameters, and is comparatively simple to use. In this article, we present further developments in the model and evaluate its new predictions compared to observed phenomena.

Seismic events can be arbitrarily divided into two main categories: regular earthquakes (EQs) and slow slip events (SSEs). They result from shear stress relaxation accumulated on the boundaries of the earth’s faults. Other common terms for SSEs are ‘slow EQs’, ‘silent EQs’, ‘aseismic strain transients’, and ‘creep events.’ Regular EQs have obviously been known for a long time [9], but SSEs were discovered comparatively recently [10-12] and are now widely observed along major transform and subduction zones [13-19]. Although large regular EQs are the most intriguing part of fault dynamics, SSEs are viewed as increasingly important, since a

large fraction of plate motion and energy release is due to aseismic movement along a fault [20]. The various mechanisms of stress relaxation therefore need to be understood in greater detail [21].

The similarities and differences between regular EQs and SSEs provide important observational constraints on any theory describing their underlying mechanisms. Both release substantial amounts of accumulated strain energy by spontaneously initiated slip along a fault. However, the latter does not radiate seismic waves. Whereas the rupture propagation velocity of EQs is about 3 km/s3 \mathrm{~km} / \mathrm{s} and the slip velocity is about 1 m/s1 \mathrm{~m} / \mathrm{s}, the propagation speed of SSEs is less by up to five orders of magnitude and the slip velocity is less by up to ten orders of magnitude. The scaling laws are also different [22-24].

The focal area of an EQ/SSE is usually located at depths inaccessible to direct observation. Seismic waves and/or transient surface deformation, observed on the Earth’s surface or at shallow depth, are the only manifestations of seismic events. Consequently, modeling is especially important for understanding processes at their source. However, the phenomenology underlying seismic events is sufficiently complicated that there are no governing equations that adequately describe crust rheology, composition, and structure, as well as the highly nonlinear character of frictional sliding along faults [25]. Therefore, various approaches have been developed to construct models for specific purposes.

The early models had exact analytical solutions that were applied to the dynamics of a crack in elastic solids. This approach revealed the fundamental relations between kinematic variables (rupture velocity, slip velocity, displacement) and dynamic variables (shear and normal stress, stress drop) characterizing EQs [26-30].

The limitations of purely analytical models were subsequently elucidated [28], leading to the development of computational models that solve a system of elasto-dynamic equations with boundary conditions given as a frictional law on crack boundaries, e.g. [31-37]. Whereas these essentially 3D models are able to describe some basic characteristics of an earthquake and accompanying seismic radiation, they are beset by various challenges, both mathematical (singularities at the tip of the crack; accuracy and relevance of the chosen boundary conditions, i.e., fault constitutive relations) and computational (issues of convergence, calculational load, and efficiency). Concurrently, a simpler 1D mass-spring model was presented [38]. This Burridge-Knopoff model has been used primarily to describe the collective behavior and

statistical features of EQs. It reproduces the major empirical laws of observed seismicity, i.e., earthquake recurrence, the Gutenberg-Richter law, foreshock and aftershock activities, e.g. [3944]. In general, however, the Burridge-Knopoff model is not adequate to describe the dynamics of an individual event.

A different methodology developed by Tse and Rice [45] is based on what is known as the rate-state friction law discovered in laboratory experiments by Dieterich [46] and Ruina [47]. The main advantages of this widely used approach are (1) its simplicity (the system of two basic equations describing the rate-state frictional law can be solved analytically in the 1D approximation and is comparatively easy to analyze computationally in the 2D approximation) and (2) its universality (the model can be used to describe both EQs, e.g. [48-51], and SSEs [5258]). However, the model can only provide a simplistic picture of rupture dynamics, which are solely determined by the ratio of the rate-state constants aa and bb. Moreover, the model predicts that seismic events occur only for b>ab>a, which is the regime of velocity-weakening frictional sliding. Yet, there are indications, both experimental [59] and theoretical [8, 60, 61], that frictional sliding instabilities are not confined only to velocity-weakening rocks.

We developed a model that is also suitable for the description of both EQs and SSEs [1-8]. The advantages of this model are (1) it is an intrinsically nonlinear,dynamical model, rooted in the Newtonian equations of motions; (2) all model parameters have explicit and unambiguous physical correlates; (3) it describes processes over a wide range of conditions, from very fast processes such as regular EQs down to very slow processes such as silent EQs with a few month rise time; (4) there are analytical solutions to the model for certain initial conditions; (5) the model is not computationally intensive.

The first application of our model was to a general description of the dynamics of crustal faults [1]. This approach provided quantitative values for the dynamic variables describing plate boundary EQs as well as slow earthquakes and afterslip and explained why temporal and spatial distributions of aftershocks obey, respectively, the Omori and a 1/r1 / \mathrm{r} laws. The model was then applied to ‘episodic tremor and slip’ (ETS) phenomenon and successfully described the basic morphological features of ETSs [2, 4, 5]. The next application was devoted to the description of the transition process from static to dynamic frictional regimes [3]. The significant result of that article was the demonstration that the rupture velocity is defined by the shear to normal stress ratio, in quantitative argument with laboratory experiments [62]. In subsequent articles [6, 7], the

governing equations were modified to include a simple velocity-dependent frictional term. Counterintuitively, the model showed that the rupture velocity does not depend on friction. Recently, we re-examined the necessary and sufficient conditions for instability of frictional sliding in the framework of our model, incorporating the rate-and-state frictional law [8]. The analysis used the standard technique of linearizing a set of equations derived from the model and investigating their behavior under small perturbations. We showed that sliding instability can develop in velocity-strengthening as well as in velocity-weakening frictional sliding, in contrast to the prevailing view that only velocity-weakening sliding leads to instability. The onset of instability is a separate topic, however. Here, we are interested in seismic dynamics assuming an instability has already developed, i.e., after nucleation of the event.

The model, to date, enables one to simulate spatial-temporal distributions of stress, slip, slip velocity, and rupture propagation velocity. There are multiple evidences, both experimental [6264] and theoretical [65-70], suggesting that stress heterogeneity plays a crucial role in nucleation and dynamics of the rupture. The spatial distribution and magnitude of heterogeneity are the main factors that determine the dynamics of a seismic event. A distinguishing feature of our model is that the stress heterogeneity is introduced as an initial condition. We have used only a simple step-like distribution previously.

In the present study, we consider three different initial conditions that represent realistic stress heterogeneities for transform and subduction faults. Moreover, rather than using our previously simple velocity-dependent frictional term, our model now incorporates the rate-state [46,47][46,47] and flash heating [71-73] frictional laws. We obtain the conditions that determine whether a seismic event becomes an EQ or a SSE by calculating the dependence of rupture characteristics (propagation velocity, slip velocity, displacement, and stress drop) on the stress spatial distribution and the coefficients in the frictional law used. In addition, the discovery of self-healing slip pulses [74] motivates a search for the mechanisms and/or conditions generating a pulse-like mode rather than a crack-like mode [44, 56, 75-86]. We determine such conditions in the framework of our model.

In this paper, we answer the following questions:

• What is the difference between EQs and SSEs from the point-of-view of their physical mechanisms?

• What are the critical parameters and their values which determine the precise spatial and temporal distribution of rupture velocity, slip velocity, slip, and stress drop?

• What is the role, both qualitative and quantitative, of stress heterogeneity in EQ and SSE development?

• What conditions (quantitative) determine the appearance of pulse-like ruptures versus crack-like ruptures during EQs and SSEs?

• What is the role of the rate-state frictional law in the evolution of seismic events?

Section II describes the model and the system of governing equations. Section III presents the results of simulations for the three initial conditions considered. The results are discussed in the context of known conventional results. The final section summarizes our findings.

II. Model

Our approach to modeling EQs and SSEs derives from the Frenkel-Kontorova (FK) model [87], as described in our previous publications [1, 3, 6, 7, 88]. Here, we briefly recap the derivation and introduce extensions to the model that are the topic of this article. We first describe the FK model that is used as the foundation. This is followed in section B by modifications we made to the FK model to adapt it to the study of frictional sliding. The sineGordon (sG) equation which naturally emerges is relevant to the description of slip dynamics at such a small scale as to be not only unobservable, but unnecessary for understanding the macroscopic phenomena of interest. Section C therefore outlines our derivation of sG modulation equations that are more suitable and more relevant to the description of large-scale phenomena. Sections D and E introduce the new elements incorporated in the current model. Table 1 provides values of the parameters used in the simulations.

A. Frenkel-Kontorova model

The FK model was proposed to explain why plasticity in crystalline materials is initiated by stress that is a few orders of magnitude smaller than was expected from prevailing analyses. It assumes the presence of defects in the crystal lattice. The existence of such defects, linear dislocations, was indeed confirmed later [89]. The reason that dislocations dramatically reduce the shear stress threshold for the onset of plasticity is that the movement of a defect requires

much less shear stress than the stress necessary to overcome the potential energy barrier between an atom and the lattice. The movement of a dislocation is accompanied by a local relative shift of the crystal parts. Thus, the plastic deformation occurs due to movement of multiple dislocations. The analogy would be to caterpillar-like motion. This process may be described by a nonlinear equation [89]:
M∂2ui∂t2−Kb(ui+1−2ui+ui−1)+Fdsin⁡2πb0ui=Fi−fiM \frac{\partial^{2} u_{i}}{\partial t^{2}}-K_{b}\left(u_{i+1}-2 u_{i}+u_{i-1}\right)+F_{d} \sin \frac{2 \pi}{b_{0}} u_{i}=F_{i}-f_{i},
where uiu_{i} is the shift of atom ii relative to its equilibrium position in the direction, xx, of plastic deformation, tt is time, b0b_{0} is a distance between atoms, MM is the mass of the atom, FdF_{d} is the amplitude of the periodic force occurring due to the presence of atomic layers beneath and above the plane of plastic deformation, FiF_{i} is the external force, and fif_{i} is the frictional (dissipative) force. The atoms interact with the nearest neighbors on either side via harmonic spring forces of constant KbK_{b}. Assuming that the crystal material has a volume density ρ\rho and an interatomic distance b0b_{0} in all directions, we can express the coefficients in Eq. (1) as (see [89] for details): M=ρb03,Kb=2Gb01−v,Fd=Gb022πM=\rho b_{0}^{3}, K_{b}=\frac{2 G b_{0}}{1-v}, F_{d}=\frac{G b_{0}^{2}}{2 \pi}, where GG is the shear modulus and vv is the Poisson ratio. Equation (1) can then be presented in the form
∂2(2πu/b0)∂(tc/b0)2−∂2(2πu/b0)∂(x/b0)2+1−v2sin⁡2πb0u=(σx−σf)2π(1−v)2G\frac{\partial^{2}\left(2 \pi u / b_{0}\right)}{\partial\left(t c / b_{0}\right)^{2}}-\frac{\partial^{2}\left(2 \pi u / b_{0}\right)}{\partial\left(x / b_{0}\right)^{2}}+\frac{1-v}{2} \sin \frac{2 \pi}{b_{0}} u=\left(\sigma_{x}-\sigma_{f}\right) \frac{2 \pi(1-v)}{2 G},
where the second term on the left hand side is obtained in the continuum limit approximation, σx\sigma_{x} is shear stress, σf\sigma_{f} is frictional force per unit area, and c2=2Gρ(1−v)≡ci2(1−2v)(1−v)2c^{2}=\frac{2 G}{\rho(1-v)} \equiv \frac{c_{i}^{2}(1-2 v)}{(1-v)^{2}}, given in terms of the longitudinal acoustic velocity, cic_{i}. Note that cc is bounded by clc_{l} and the shear wave velocity, csc_{s}.

B. FK modified for frictional slip dynamics

We apply the FK model to describe slip dynamics by modeling macroscopic friction as the result of interactions between asperities on the surfaces of two sliding solids. We consider the asperities on one of the frictional surfaces as forming a linear chain of equally spaced balls. The

simplest model assumes the asperities on the opposite frictional surface form a rigid fixed substrate with the spacing between asperities equal to the ball spacing. The variable uu then represents the displacement of a given ball at position xx and time tt relative to its non-displaced position defined as being in the middle between two of the opposite asperities.

The maximum force between surfaces in the case of plasticity, FdF_{d}, becomes (1−v)/2(1-v) / 2 in Eq. 2. For modeling frictional sliding, the value of this coefficient should be reduced by the ratio between nominal and actual contact areas, which has been shown [90, 91] to be equal to the ratio ΣN/σp\Sigma_{N} / \sigma_{p}, where ΣN\Sigma_{N} is the effective normal stress (i.e., normal stress minus pore pressure) and σp\sigma_{p} is the penetration hardness. The hardness of the polycrystalline materials and shear modulus are related by a simple formula, σp≈G/(2π)\sigma_{p} \approx G /(2 \pi) [92], hence the ratio ΣNσp≡2πΣNG\frac{\Sigma_{N}}{\sigma_{p}} \equiv \frac{2 \pi \Sigma_{N}}{G}. We also have σf=μΣN\sigma_{f}=\mu \Sigma_{N} in terms of the standard coefficient of friction, μ\mu. Then, to describe frictional sliding, Eq. (2) becomes
∂2(2πu/b0)∂(tc/b0)2−∂2(2πu/b0)∂(x/b0)2+2πΣNG1−v2sin⁡2πb0u=(σs−σf)2π(1−v)2G\frac{\partial^{2}\left(2 \pi u / b_{0}\right)}{\partial\left(t c / b_{0}\right)^{2}}-\frac{\partial^{2}\left(2 \pi u / b_{0}\right)}{\partial\left(x / b_{0}\right)^{2}}+\frac{2 \pi \Sigma_{N}}{G} \frac{1-v}{2} \sin \frac{2 \pi}{b_{0}} u=\left(\sigma_{s}-\sigma_{f}\right) \frac{2 \pi(1-v)}{2 G},
where now b0b_{0} is the typical distance between asperities. Dividing both sides by 2πΣNG1−v2\frac{2 \pi \Sigma_{N}}{G} \frac{1-v}{2} and introducing dimensionless variables u~\tilde{u} in units of b0/(2π),x~b_{0} /(2 \pi), \tilde{x} in units of [2G2πΣN(1−v)]1/2b0,t~\left[\frac{2 G}{2 \pi \Sigma_{N}(1-v)}\right]^{1 / 2} b_{0}, \tilde{t} in units of [2G2πΣN(1−v)]1/2b0c\left[\frac{2 G}{2 \pi \Sigma_{N}(1-v)}\right]^{1 / 2} \frac{b_{0}}{c}, and σ~x~\tilde{\sigma}_{\tilde{x}} and σ~f=μ\tilde{\sigma}_{f}=\mu in units of ΣN\Sigma_{N}, we obtain the well-known sineGordon (sG) equation:
∂2u~∂t~2−∂2u~∂x~2+sin⁡u~=σ~s−σ~f\frac{\partial^{2} \tilde{u}}{\partial \tilde{t}^{2}}-\frac{\partial^{2} \tilde{u}}{\partial \tilde{x}^{2}}+\sin \tilde{u}=\tilde{\sigma}_{s}-\tilde{\sigma}_{f},

C. The sine-Gordon modulation equation

There are three fundamental solutions of the sG equation: kinks (or solitons), plasma waves, and breathers. In the framework of plasticity, a soliton represents a linear dislocation. In the context of friction, we model a soliton as an area of accumulated, enhanced stress nucleated on the surface by shear stress applied to an inhomogeneous surface. We will use the term

“macroscopic dislocation” to designate such areas. In crystals, the displacement of a dislocation under the applied shear stress (its “mobility”) is much larger than the corresponding translation of the entire atomic plane. In the same way, the mobility of a macroscopic dislocation over the “bumpy road” of a frictional surface (i.e., the quasi-periodic substrate) is larger than the mobility of the whole surface. In both cases, the displacement of a dislocation (a pre-stressed area) requires much less external stress. Therefore, in our model, the relative sliding of two bodies along a planar interface, which occurs due to the movement of macroscopic dislocations, can be initiated by external stress smaller than μΣS\mu \Sigma_{S}, which is consistent with laboratory experiments [62]. We note that a dislocation may propagate with any velocity ranging from zero to cc.

The sG equation has been widely used to describe nanoscale friction (see, e.g., [93] and references therein). It has also been applied to macroscopic friction [94-96], but it is not particularly useful for this case. The sG equation is relevant for quantitatively describing the dynamics of only a few dislocations. It is not useful for describing plastic deformation in a macroscopic crystal, which may include millions of dislocations per centimeter. The same is true for macroscopic friction. If one requires a detailed description of the relative movement of two bodies at a scale of a few asperities, the sG equation is the right tool. However, for objects which include millions of asperities per cm2\mathrm{cm}^{2} (i.e., thousands of macroscopic dislocations), the pure sG is not an appropriate instrument. The details it provides of dynamics at the microscopic level are not directly observable in standard macroscopic experiments. The actual position of any particular dislocation is not important.

The dynamics of a sequence of a large number of dislocations are better described by variables averaged in space and time over the dislocation periods. Witham developed a technique for transforming an equation describing the behavior of given variables to a system of equations describing the behavior of the variables averaged over space and time [97]. Applying this method to the sG equation, we derived a system of sG modulation equations [88], expressed again in dimensionless variables, which is suitable for describing macroscopic friction:
−4πK[ m(1−U2)]1/2[UtηU2−1+mtU2m+UxUηU2−1+mx12m]=ΣS−μ-\frac{4}{\pi} \frac{K}{\left[\mathrm{~m}\left(1-U^{2}\right)\right]^{1 / 2}}\left[U_{t} \frac{\eta}{U^{2}-1}+m_{t} \frac{U}{2 m}+U_{x} \frac{U \eta}{U^{2}-1}+m_{x} \frac{1}{2 m}\right]=\Sigma_{S}-\mu
UtUU2−1+mtη2mmt+Ux1U2−1+mxUη2mmt=0U_{t} \frac{U}{U^{2}-1}+m_{t} \frac{\eta}{2 m m_{t}}+U_{x} \frac{1}{U^{2}-1}+m_{x} \frac{U \eta}{2 m m_{t}}=0

Here mm is the modulus of the Jacoby elliptic function (0<m≤1),m1=1−m,η=E/K(0<m \leq 1), m_{1}=1-m, \eta=E / K, where K(m)K(m)K(m) K(m) and E(m)E(m)E(m) E(m) are the complete elliptic integrals of the first and second kind, respectively, UU is the wave velocity in units of velocity cc, with −1≤U≤1-1 \leq U \leq 1. Subscripts xx and tt denote the derivative with respect to the associated variable. The friction coefficient μ\mu and (dimensionless) shear stress Σs\Sigma_{s} now represent the more relevant macroscopic averages over the dislocation periods.

D. Relation to macroscopic physical variables

The variables m(x,t)m(x, t) and U(x,t)U(x, t) obtained by solving Eqs. (5) can be associated with macroscopic physical variables of interest. The (dimensionless) shear stress and slip velocity WW (in units of c(ΣN2πG1−v2)1/2c\left(\frac{\Sigma_{N}}{2 \pi G} \frac{1-v}{2}\right)^{1 / 2} are given by [1]
ΣS=πK[m(1−U2)]1/2\Sigma_{S}=\frac{\pi}{K\left[m\left(1-U^{2}\right)\right]^{1 / 2}},
W=UΣSW=U \Sigma_{S}.
Slip SS, which is directly related to the suitable average for uˉ(x,t)\bar{u}(x, t) in Eq. (4) that motivated the derivation of Eqs. (5), is simply
S=∫0tWdtS=\int_{0}^{t} W d t,
while stress drop, ΔΣ\Delta \Sigma, is
ΔΣ(x,t)=ΣS(x,t)−ΣS(x,0)\Delta \Sigma(x, t)=\Sigma_{S}(x, t)-\Sigma_{S}(x, 0).
Given initial values at t=0t=0 for the spatial distribution of ΣS\Sigma_{S} and WW, one can use Eqs. (6) and (7) to obtain initial values for m(x,0)m(x, 0) and U(x,0)U(x, 0). One must also input an initial value for friction, which is the topic of the next section. Equations (5) can then be written using only the two variables mm and UU and solved numerically. Equations (6) and (7) then give the redistribution of shear stress, ΣS(x,t)\Sigma_{S}(x, t) and slip velocity, W(x,t)W(x, t), from which are derived the slip and stress drop using Eqs. (8) and (9). Any of these spatial and temporal profiles will show the evolving rupture front or edge as a clear discontinuity delineating the modified values in the rupture at

times tt from the initial values at positions ahead of the front. Then the average rupture velocity, VV, is simply the distance the rupture front has travelled in the corresponding time interval.

E. Friction

In previous work, our model has either not included friction or simply set it proportional to slip velocity [6, 7]. In the present work, we incorporate a more physically realistic friction into our model. For slip velocity W≤0.1 m/sW \leq 0.1 \mathrm{~m} / \mathrm{s}, we will use the well-known Dieterich-Ruina rate-state model [46,47][46,47] for the coefficient of friction, μ\mu :
μ(W,θ)=μ0+aln⁡(WW0)+bln⁡(W0θDc)\mu(W, \theta)=\mu_{0}+a \ln \left(\frac{W}{W_{0}}\right)+b \ln \left(\frac{W_{0} \theta}{D_{c}}\right)
dθdt=1−WθDc\frac{d \theta}{d t}=1-\frac{W \theta}{D_{c}},
where aa and bb are empirical dimensionless constants, θ\theta is the state variable in units of b0G/(2πΣNc)b_{0} G /\left(2 \pi \Sigma_{N} c\right), and DcD_{c} is the characteristic slip distance. The connection between parameters b0b_{0} and DcD_{c} is unknown. We assume that b0/(2πDc)≈1b_{0} /\left(2 \pi D_{c}\right) \approx 1. The latter is true if asperities have a sinusoidal shape and if DcD_{c} is equal to the root-mean-square roughness of the surface. When θ˙=\dot{\theta}= 0 , we have the steady state values θSS=Dc/W\theta_{S S}=D_{c} / W, in which case, the steady state friction μSS\mu_{S S} is μ(W,θss)≡μss(W)=μ0+(a−b)ln⁡(WW0)\mu\left(W, \theta_{s s}\right) \equiv \mu_{s s}(W)=\mu_{0}+(a-b) \ln \left(\frac{W}{W_{0}}\right),
giving μSS(W0)=μ0\mu_{S S}\left(W_{0}\right)=\mu_{0}. Representative reference values for μ0\mu_{0} and W0W_{0} are provided in Table 1.
For W>0.1 m/sW>0.1 \mathrm{~m} / \mathrm{s}, flash heating leads to a low frictional strength [69-71]. To describe friction for slip velocities W>Wflash =0.1 m/sW>W_{\text {flash }}=0.1 \mathrm{~m} / \mathrm{s}, Rice [98] proposed a simple formula:
μ(W)=[μss(W=Wflash )−μflash ]Wflash W+μflash \mu(W)=\left[\mu_{s s}\left(W=W_{\text {flash }}\right)-\mu_{\text {flash }}\right] \frac{W_{\text {flash }}}{W}+\mu_{\text {flash }},
where μflash \mu_{\text {flash }} is the minimum value of the friction coefficient resulting from flash heating, determined from experiments to be ∼0.2\sim 0.2 [71-73]. However, the effective friction also tends to increase due to radiation damping (related to the power of the radiated seismic waves) as the velocity WW increases. To include this competing effect, we modify the previous relation to
μ(W)=[(μss(Wflash )−μflash )Wflash W+μflash ][1+α(W−Wflash Wflash )2]\mu(W)=\left[\left(\mu_{s s}\left(W_{\text {flash }}\right)-\mu_{\text {flash }}\right) \frac{W_{\text {flash }}}{W}+\mu_{\text {flash }}\right]\left[1+\alpha\left(\frac{W-W_{\text {flash }}}{W_{\text {flash }}}\right)^{2}\right],

where α=(Wflash /Wm)2(μ0−μflash )(1−Wflash /Wm)[(μ0−μflash )Wflash /Wm+μflash ]\alpha=\frac{\left(W_{\text {flash }} / W_{\mathrm{m}}\right)^{2}\left(\mu_{0}-\mu_{\text {flash }}\right)}{\left(1-W_{\text {flash }} / W_{\mathrm{m}}\right)\left[\left(\mu_{0}-\mu_{\text {flash }}\right) W_{\text {flash }} / W_{\mathrm{m}}+\mu_{\text {flash }}\right]}, and (dimensionless) Wm=1.1×10−2W_{m}=1.1 \times 10^{-2} (equal to 1 m/s1 \mathrm{~m} / \mathrm{s} in physical units) is a representative average value for the maximum slip velocity within a rupture for an EQ event. This expression for α\alpha is derived to satisfy the condition μ(Wm)=μss(Wflash )\mu\left(W_{m}\right)=\mu_{s s}\left(W_{\text {flash }}\right). An example of the dependence of friction on slip velocity, μ(W)\mu(W), is shown in Fig. 1 for the particular value θ=θSX\theta=\theta_{S X} (steady state). For the velocity weakening case chosen ( a<ba<b, see Table 1), friction steadily decreases, as given by Eqs. (10), for slip velocities W≤Wflash W \leq W_{\text {flash }}. Then as WW increases to W>Wflash W>W_{\text {flash }}, friction abruptly decreases further and finally sharply increases with further increase of WW, as given by Eq. (11).

F. Initial conditions used in the model

To solve the system of Eqs. (5) - (11), two initial conditions must be specified. The initial slip velocity is chosen to be the average slip velocity along a fault, W0=W(x,t=0)W_{0}=W(x, t=0). A typical value is W0=3 cm/yrW_{0}=3 \mathrm{~cm} / \mathrm{yr}. In contrast to conventional approaches, we will use spatially heterogeneous shear stress as an initial condition. The reason for this is that sliding instabilities develop in locations with a gradient in the ratio of shear to effective normal stress [8]. Hereafter we will use the term “stress ratio” instead of “shear to effective normal stress ratio”.

Three configurations of the initial stress distribution will be used for the simulations (shear stress applied, without loss of generality, in the positive xx-direction):
ΣS(x,t=0)=ΣS0−δΣπtan⁡−1(βx)\Sigma_{S}(x, t=0)=\Sigma_{S 0}-\frac{\delta \Sigma}{\pi} \tan ^{-1}(\beta x)
ΣS(x,t=0)=ΣS0−δΣpπtan⁡−1(βx)e−ζx\Sigma_{S}(x, t=0)=\Sigma_{S 0}-\frac{\delta \Sigma_{p}}{\pi} \tan ^{-1}(\beta x) e^{-\zeta x}
ΣS(x,t=0)=ΣS0+δΣ2cosh⁡(ζx)\Sigma_{S}(x, t=0)=\Sigma_{S 0}+\frac{\delta \Sigma}{2 \cosh (\zeta x)},
where ΣS0\Sigma_{S 0} is the value of the average shear stress in the distribution, δΣ\delta \Sigma is the difference between the maximum and minimum stress, β\beta regulates the slope of the profile connecting the maximum and minimum stress values, and ζ\zeta determines the spatial extent of variations in the stress distribution. In order to properly compare the simulation results, the maximum stress in all

three cases is set equal, which requires Eb=π2e−∣xmax⁡∣tan⁡−1(β∣xmax⁡∣)E\mathfrak{E}_{b}=\frac{\pi}{2} \frac{e^{-\left|x_{\max }\right|}}{\tan ^{-1}\left(\beta\left|x_{\max }\right|\right)} \mathfrak{E} in Eq. (12b), where xmax⁡x_{\max } is the location of the maximum in the stress profile.

Spatial profiles for the three initial conditions are illustrated in Figs. 2(a-c). The first two configurations mimic the spatial distribution of shear stress around an obstacle, which prevents relative sliding of the plates along a fault. Stress accumulates on one side of the obstacle and there is a stress deficit on the other side. The third configuration approximates effects due to localized reductions in the effective normal stress, due, for example, to changes in pore pressure or “nonplanar interfaces being compressed into full contact” (e.g. [99]).

A FORTRAN code was developed to obtain numerical solutions. The code allows one to study the evolution of frictional processes employing the predictor-corrector scheme. Typical values for the effective normal stress used in the model are 60 MPa for regular EQs and 100 MPa for SSEs, based on the lithostatic pressure at depth hh, equal to ρgh\rho g h.

III. Results and Discussion

New results from the updated model are presented here. In Section A, we show that the model is capable of simulating a wide range of seismic events, from a typical regular EQ with slip velocity ∼1 m/s\sim 1 \mathrm{~m} / \mathrm{s} and rupture velocity ∼cs\sim c_{s} to a SSE with slip velocity ∼1 m/\sim 1 \mathrm{~m} / year and rupture velocity ∼1 km/\sim 1 \mathrm{~km} / day. Then we determine the specific conditions (initial stress spatial distribution and frictional parameters) that produce either an EQ or SSE. In Section B, we provide examples to show the model can account for both crack-like and pulse-like ruptures. We then identify the conditions that control the rupture mode. In Section C, we consider additional differences in rupture shape that result from differences in the initial stress heterogeneity.

A. Conditions for regular EQ versus SSE

Let us first consider a particularly simple example of stress heterogeneity in which the initial stress is the non-localized step-like function given by Eq. 12(a) and depicted in Fig. 2(a). We assume that a rupture (seismic event) is initiated at t=0t=0 and at the position x=0x=0 where the gradient in the stress ratio is maximum.

Values ΣS0=0.7\Sigma_{S 0}=0.7 and δΣ=0.1\delta \Sigma=0.1 give results that are typical for a regular EQ. Figure 3 shows the resulting solutions for the spatial distribution of shear stress, ΣS\Sigma_{S}, and slip velocity, WW, at three different times. Further details of the spatial-temporal distribution of ΣS\Sigma_{S} and WW, as well as stress drop, δΣ\delta \Sigma, and slip, SS, are shown in Fig. 4. Recall that the average rupture velocity, VV, is the distance the front travels (the location of the rupture edge, as shown in the figures) divided by the time travelled. As one can see, the rupture propagates in two directions with slightly different velocities that we calculate as V−=3.30 km/s=0.9csV^{-}=3.30 \mathrm{~km} / \mathrm{s}=0.9 c_{\mathrm{s}} to the left and V+=2.24 km/s=0.61csV^{+}=2.24 \mathrm{~km} / \mathrm{s}=0.61 c_{\mathrm{s}} to the right. The rupture propagates faster in the region of greater shear stress, V−>V+V^{-}>V^{+}, consistent with our previous finding [1, 3]. As can be seen from Fig. 4, WW quickly reaches a value of ∼1 m/s\sim 1 \mathrm{~m} / \mathrm{s} near the initiation point x=0x=0 and remains constant at this point thereafter. In the rest of the rupture area, WW is ∼0.8 m/s\sim 0.8 \mathrm{~m} / \mathrm{s}. The slip is maximum at the initiation point and reaches ∼10 m\sim 10 \mathrm{~m} after ∼10sec\sim 10 \mathrm{sec}. The maximum stress drop shown in the figure is about 3 MPa , which is about 7%7 \% of the average initial shear stress of ∼42MPa\sim 42 \mathrm{MPa}.

Performing the same calculations using the values ΣS0=0.18\Sigma_{S 0}=0.18 and δΣ=5×10−4\delta \Sigma=5 \times 10^{-4} produces a SSE. The results, converted to physical units using the SSE normalization factor ΣN=100MPa\Sigma_{N}=100 \mathrm{MPa}, are shown in Fig. 5. The slip velocity W≈2.3 m/W \approx 2.3 \mathrm{~m} / year, and the rupture velocities are V−≈V+≈0.88 km/V^{-} \approx V^{+} \approx 0.88 \mathrm{~km} / day, which are, respectively, seven and five orders of magnitude less than slip and rupture velocities of a regular EQ. The stress drop is about 40 KPa , two orders of magnitude less than for a regular EQ. The stress drop could be even smaller than given by our example using ΣN=100MPa\Sigma_{N}=100 \mathrm{MPa}, since there are indications [4,5][4,5] that the effective normal stress in subduction zones may be much less than the lithostatic pressure.

Thus, the differences between the examples illustrated in Figs. 4 and 5 indicate that the dynamic variables describing the rupture depend critically on the values of ΣS0\Sigma_{S 0} and δΣ\delta \Sigma relative to ΣN\Sigma_{N}. To examine this dependence in detail, we solved Eqs. (5−12)(5-12) for wide ranges of ΣS0\Sigma_{S 0} (from 0.15 to 0.9 ) and δΣ\delta \Sigma (from 0.001 to 0.1 ). The initial stress distribution of Eq. 12a, Fig. 2(a) and friction, Eqs. (10) and (11), were inputs to the model. They were calculated using parameters listed in Table 1, which gives the velocity-weakening (a<b)(a<b) profile for friction (see Fig. 1). Results are plotted in Fig. 6. The maximum calculated slip velocity in the rupture, Wmax⁡W_{\max }, is plotted in the upper panel for different values of δΣ\delta \Sigma as a function of the maximum stress

Σmax⁡=Σ50+δΣ/2\Sigma_{\max }=\Sigma_{50}+\delta \Sigma / 2, normalized to ΣN\Sigma_{N}. The solid curves are the results of the friction model and show the effects of changing from rate-state, Eq. (10) to flash heating, Eq. (11), for values of Σmax⁡\Sigma_{\max } that give Wmax⁡≥Wflash W_{\max } \geq W_{\text {flash }}. The dashed curves are the results of using only rate-state friction for all Σmax⁡\Sigma_{\max }. The rupture velocity V−V^{-}V−of the front propagating in the direction of decreasing xx is plotted in the lower panel.

Changing Σmax⁡/ΣN\Sigma_{\max } / \Sigma_{N} by less than an order of magnitude changes Wmax⁡W_{\max } by about nine orders of magnitude and V−V^{-}V−by about seven orders. When Wmax⁡=Wflash W_{\max }=W_{\text {flash }}, friction abruptly decreases due to flash heating, and the (dimensionless) slip velocity jumps to values Wmax⁡∼10−2W_{\max } \sim 10^{-2}, (see abrupt steps in the solid curves of the upper panel) corresponding to physical values ≥1 m/s\geq 1 \mathrm{~m} / \mathrm{s} that are characteristic of EQs. There is also a jump in V−forthesametransitionvaluesofforthesametransitionvaluesofSigmamax/SigmaNV^{-}forthesametransitionvaluesoffor the same transition values of \Sigma_{\max } / \Sigma_{N}V−forthesametransitionvaluesofforthesametransitionvaluesofSigmamax/SigmaN, as shown in the lower panel. If flash heating is not considered and only rate-state friction is used for all WW, the slip and rupture velocities are much smaller with increasing Σmax⁡/ΣN\Sigma_{\max } / \Sigma_{N} (see dotted curves in the upper panel and curves in the bottom right rectangle of the lower panel). Almost identical results (not shown) are obtained for the other two initial stress profiles, Eqs. 12b12 b and cc. These results do not appear to depend sensitively on the exact profile of the shear stress.

Thus, we conclude that in our model, the type of seismic event, i.e. EQ or SSE, is determined primarily by the stress ratio Σmax⁡/ΣN\Sigma_{\max } / \Sigma_{N}. Furthermore, EQs, i.e., slip with maximum slip velocity in the rupture Wmax⁡∼1 m/s(∼0.01W_{\max } \sim 1 \mathrm{~m} / \mathrm{s}(\sim 0.01 in dimensionless units), only occur if there is a dramatic jump in slip velocity due to a pronounced decrease in friction due to flash heating or other mechanisms when the velocity reaches Wflash W_{\text {flash }}. Rate-state friction alone does not provide the necessary slip velocity using experimental values for the parameters aa and bb unless unrealistically large values Σmax⁡≥0.9ΣN\Sigma_{\max } \geq 0.9 \Sigma_{N} are used.

In previous work, we showed that, in contrast to common assumption, sliding instability can develop for velocity-strengthening ( a>ba>b ), as well as velocity-weakening ( a<ba<b ), frictional sliding [8]. Instability in velocity-strengthening sliding was found to be the result of a sufficiently large gradient in the stress ratio. Velocity-weakening instability is, nonetheless, more probable since nucleation lengths are smaller than for velocity-strengthening materials. Here, for W<Wflash W<W_{\text {flash }}, we find in addition that when the rupture is already nucleated, i.e., an instability has developed, the dynamics of the rupture are the same for the velocity-weakening and velocity-

strengthening cases. In other words, the dynamics of the rupture do not depend on the relative values of aa and bb in the rate-state law. The foregoing holds for all SSEs since W<Wfash W<W_{\text {fash }} in all these seismic events.

B. Conditions for crack-like versus self-healing pulse rupture

In the previous section, we found that the primary determinants of whether a rupture event manifests as an EQ or a SSE are the following values used as inputs to the model for the initial stress distribution: the maximum value in shear to normal stress ratio, Σmax⁡/ΣN\Sigma_{\max } / \Sigma_{N}, and δΣ\delta \Sigma, the difference between the maximum and minimum stress. The results were not sensitive to other details of the stress profile, so we focused on examples obtained using the initial stress distribution of Eq. 12a. Having established that the model can account for both EQs and SSEs, as well as providing the conditions that distinguish between them (e.g., Fig. 6), we now look for rupture properties that do depend on particulars of the spatial heterogeneity in the initial stress profile.

Note that all rupture examples considered so far have been crack-like, i.e., the slip velocity at any point within the rupture is approximately constant in time. In the context of a seismic event it means that the rise time equals the slip time at a rupture initiation point. Heaton [74] showed that for some EQs the slip time at some point of a rupture is much less than the EQ rise time indicating that at least in these particular areas the rupture takes the form of a self-healing pulse.

Consider the spatially localized inhomogeneity in the initial stress given by Eq. 12b12 b (see also Fig. 2b2 b ). The results of simulations using values (Σ50=0.7,δΣ=0.17)\left(\Sigma_{50}=0.7, \delta \Sigma=0.17\right) that give an EQ are shown in Figs. 7 and 8. We find a crack-like rupture that propagates at velocities VpmV^{ \pm}Vpmin opposite directions. The localized slip to the right, i.e. in the shear stress direction, causes a stress drop at the left side of the rupture and stress rise at the right side (Fig. 7 illustrates the development of this structure). As a result, the disturbance area expands in both directions. The maximum slip velocity in the rupture is about 1 m/s1 \mathrm{~m} / \mathrm{s}, as expected for an EQ, and the maximum stress drop is about 4 MPa .

The results are similar to those obtained using the initial condition of Eq. 12a (see Figs. 3 and 4) both qualitatively and quantitatively, i.e. the rupture is an expanding crack. Investigating further, we find the rupture is always crack-like for the initial stress profile of Eq. 12a, which is asymmetrically homogeneous outside the heterogeneous region. For Eq. 12b12 b, the rupture is

crack-like only if Σ50\Sigma_{50} and δΣ\delta \Sigma are large enough that the system evolves to the abrupt velocityweakening (flash heating) conditions, so that the maximum slip velocity in the rupture becomes larger than Wflash =0.1 m/sW_{\text {flash }}=0.1 \mathrm{~m} / \mathrm{s}. For smaller values of Σ50\Sigma_{50} and/or δΣ\delta \Sigma the system can behave very differently.

We consider the previous example, but reduce δΣ\delta \Sigma by a factor of 10 to produce an SSE, with results shown in Figs. 9 and 10. The rupture now expands in both directions by Heaton-like (selfhealing) pulses, i.e. the slip velocity at the rupture fronts is much larger than inside the rupture area. The maximum slip velocity is about 0.06 m/s0.06 \mathrm{~m} / \mathrm{s}, and the system never switches from rate-state to flash heating. The maximum slip and the maximum stress drop are centered at the initiation point, x=0x=0. The rupture propagation velocity is a factor 1.3 smaller than in the previously discussed case of larger maximum slip velocity (Fig. 8), where the rupture is crack-like.

Further analysis shows, in addition to W<Wflash W<W_{\text {flash }}, a pulse-like rupture requires a spatially localized heterogeneity in the initial stress profile. The solution should be obtained over a region that exceeds the extent of the heterogeneity by at least a factor of ∼3\sim 3 in order to distinguish the pulse-like character of the rupture. If flash heating is included in the simulations and Σ50\Sigma_{50} and δΣ\delta \Sigma are large enough to drive the system to the conditions where W>Wflash W>W_{\text {flash }}, then the rupture is crack-like. However, if flash heating is not included in the model, the rupture is always pulselike. The abrupt reduction of friction due to flash heating prevents the development of a pulselike rupture. The conditions that distinguish between the occurrence of crack-like and pulse-like ruptures were obtained from simulations and are summarized in Fig. 11.

C. Nuances in rupture characteristics due to stress heterogeneity

We have found the conditions that determine whether a seismic event is an EQ or SSE (see Fig. 6) and whether the resulting rupture is crack-like or pulse-like (Fig. 11). In particular, we found that these results were not sensitive to details of the stress profile, depending only on Σ50\Sigma_{50} and δΣ\delta \Sigma. In this section, we investigate effects due to different types of initial stress profile.

Consider the spatially localized inhomogeneity in the initial stress given by Eq. 12c (see also Fig. 2c). The results of simulations using values Σ50=0.7\Sigma_{50}=0.7 and δΣ=0.05\delta \Sigma=0.05 that give an EQ are shown in Fig. 12. As expected for these parameters, the rupture is crack-like. However, in contrast to the crack-like ruptures obtained previously, the region of maximum slip and

maximum stress drop is no longer centered at the rupture initiation point, x=0x=0, but shifts to the right.

For values (ΣS0=0.2,δΣ=5×10−4)\left(\Sigma_{S 0}=0.2, \delta \Sigma=5 \times 10^{-4}\right) that give a SSE (convert to physical units using ΣN=\Sigma_{N}= 100 MPa ), the rupture, shown in Fig. 13, is pulse-like, as expected. However, only one slip pulse is formed. There is no slip for points x<0x<0 even though there are still two stress pulses, one travelling to the left (opposite to the shear stress direction) and one to the right (in the direction of shear stress).

The differences in both crack-like and pulse-like behavior demonstrated using different initial stress profiles provide information on the processes underlying a seismic event. As noted, the initial conditions of Eqs. 12a12 a and 12b12 b mimic the spatial distribution of shear stress around an obstacle, which prevents relative sliding of the plates along a fault. The initial condition of Eq. 12c12 c reflects localized reductions in the effective normal stress

Our model also provides a quantitative description of episodic tremor and slip (ETS) [2, 4, 5], defined as the slow propagation of a single slip pulse along a subduction fault accompanied by massive bursts of nonvolcanic tremor. This periodic phenomenon has been observed in virtually all major subduction zones [16, 100-104]. In the Cascadia area, the pulse slip is 3-5 cm and pulse propagation velocity is ∼10 km/\sim 10 \mathrm{~km} / day [105]. In the current work, we have provided a possible mechanism for such pulses, which can be produced by choosing appropriate values of Σs\Sigma_{s} and δΣ\delta \Sigma, as shown in Fig. 13.

IV. Conclusion

The inherently nonlinear model presented here describes the spatial-temporal distribution of stress, slip velocity, slip, and rupture front propagation velocity resulting from input values of physical variables taken from experiment relevant to the study of earthquakes. There are no adjustable parameters in the model, and it is able to describe the dynamics of both regular EQs and SSEs, including ETSs. In our previous work, frictional effects were modeled using a simple velocitydependent term. The initial stress distribution was a simple step-like function. Here, we provided a more sophisticated treatment of friction and focused on the importance of stress heterogeneity in rupture nucleation and dynamics, as demonstrated in laboratory experiments [62-64] and suggested by various models [65-70]. To that end, (a) rate-state and flash-heating friction laws are

now coupled with the model’s system of governing equations and (b) three different heterogeneous spatial distributions for the shear to normal stress ratio have been considered as initial conditions.

The results have a fairly simple organizational structure that can be categorized according to the frictional law and the initial stress distribution. Within the context of these broad influences, our analysis shows that the type of seismic event (regular EQ or SSE), rupture mode (crack-like or pulse-like ruptures), and rupture velocity depend solely on (1) the averaged stress ratio (ratio of shear to effective normal stress); and (2) the gradient in this ratio.

In section III A, we first investigated implications of the frictional law used in the model. The abrupt decrease of friction, which occurs at slip velocities W≥Wflash =0.1 m/sW \geq W_{\text {flash }}=0.1 \mathrm{~m} / \mathrm{s} due to flash heating [71-73] or other mechanisms [106, 107], is crucial for the occurrence of EQs. Otherwise, only SSEs are possible. One caveat is that rate-state friction alone, using experimental values for its parameters aa and bb (see Table 1), could produce EQs if one allows unrealistically large values of the shear to normal stress (greater than ∼0.9\sim 0.9 ). Other models [86, 108] have also noted the importance of abrupt velocity-weakening friction at large slip velocity. We then quantified the maximum slip velocity, Wmax⁡W_{\max }, as a function of the stress and gradient to find the values that give Wmax⁡≥Wflash W_{\max } \geq W_{\text {flash }}, resulting in EQs, or Wmax⁡<Wflash W_{\max }<W_{\text {flash }}, resulting in SSEs (see Fig. 6).

In section III B, we next investigated the effect of spatial heterogeneity in the stress profile and found that, for localized heterogeneities as in Eqs. 12 (b) and ©, shown in Figs. 2 (b) and ©, EQs are always crack-like and SSEs are always pulse-like. We then calculated the conditions resulting in a particular rupture mode, as shown in Fig. 11. These results may be applicable to more general stress distributions. The stress profiles used in the simulations were either symmetric or anti-symmetric, and a range of localized stress distributions can be constructed from linear combinations of these two shapes. The anti-symmetric step-like distribution of Eq. 12(a), which extends indefinitely in both directions and is therefore not localized, produces only crack-like ruptures. Various models [44, 56, 75-86] have also considered the conditions leading to crack-like ruptures or self-healing (pulse-like) ruptures. Our model is distinctive in linking the type of rupture mode to spatial heterogeneity in the stress profile and providing quantitative values for the initial stress ratio and its gradient that produce a given rupture mode.

Upon further investigation of rupture modes, we found, in section III C, that a symmetric localized stress distribution always results in a single-pulse rupture mode for SSEs, which is

suggestive of ETSs. The heterogeneity profile also determines the position of the maximum slip for EQs, which, as shown, are always crack-like in our model. For the anti-symmetric stress profiles considered in Figs. 2(a) and (b), the location of the maximum displacement is at the rupture initiation point. For the symmetric stress profile shown in Fig. 2©, the maximum displacement is shifted in the positive xx-direction, which we have chosen as the direction of the applied shear stress.

Rupture velocity is also an important characteristic of EQs studied experimentally [109, 110] and by modeling [111-114]. We find that the rupture velocity for regular EQs ranges from a fraction of the shear wave velocity up to the supershear velocity. The range for SSEs is much wider, from 10−7cs10^{-7} c_{s} to 10−1cs10^{-1} c_{s}. The maximum slip velocity within the rupture for regular EQs is ≥\geq 1 m/s1 \mathrm{~m} / \mathrm{s}. For SSEs, slip velocities range widely from a few cm/\mathrm{cm} / year up to 0.1 m/s0.1 \mathrm{~m} / \mathrm{s}. Thus, our modeling of rupture and slip velocities are consistent with the observational data.

The importance of obtaining a value for the stress drop during EQs and SSEs is of considerable interest, e.g. [24, 115]. Using typical parameters for EQs and SSEs, our simulations show that stress drop varies in the range from 1%1 \% to 10%10 \% of the initial shear stress for regular EQs and from 0.1%0.1 \% to 10%10 \% for SSEs. Generally, due to differences in the initial shear stress between EQs and SSEs, the stress drop in SSEs is one to three orders of magnitude less than in EQs, consistent with other estimates [24].

Thus, various features of seismic events predicted by the model are consistent with basic results of conventional models. However, there are also important differences. (1) As shown previously [8], in the framework of our model, the development of instability is possible not only for velocity-weakening rocks, as given by conventional models, but also for velocitystrengthening rocks. The nucleation size depends on values of the rate-state parameters aa and bb, as usual, but also on the value of the gradient in the stress ratio. (2) New simulations presented here show that after an instability develops, rupture dynamics do not depend on the parameters aa and bb. It does not matter whether frictional sliding is velocity-weakening (a<b)(a<b) or velocitystrengthening (a>b)(a>b). Popular models based on the rate-state frictional law [45, 48-58] are not consistent with these results.

In summary, the critical factor determining the spatial and temporal distribution of the rupture velocity, slip velocity, slip, and stress drop, which characterize the dynamics of EQs and SSEs, is the initial spatial distribution of the stress ratio (which includes the gradient in the ratio)

prior to the onset of instability in frictional sliding. We have provided quantitative values for the distributions that determine a wide range of observed EQ and SSE phenomena. We also note that our model provides a quantitative description of certain features of ETS phenomenon (see Fig. 13 and the articles by Gershenzon et al. [2] and Gershenzon and Bambakidis [5]). However, others features of ETSs, such as rapid tremor propagation in the slip-parallel direction [105] and reverse tremor migration [116], can be described quantitatively only in 2D models, currently under development for presentation in a separate article (Part II).

ACKNOWLEDGMENT

We thank the reviewers for elucidating limitations in the submitted manuscript.

REFERENCES

[1] N. I. Gershenzon, V. G. Bykov, and G. Bambakidis, Strain waves, earthquakes, slow earthquakes, and afterslip in the framework of the Frenkel-Kontorova model, Phys. Rev. E 79, 056601 (2009). doi:10.1103/PhysRevE.79.056601.
[2] N. I. Gershenzon, G. Bambakidis, E. Hauser, A. Ghosh, and K. C. Creager, Episodic tremors and slip in Cascadia in the framework of the Frenkel-Kontorova model. Geophys. Res. Lett. 38, L01309 (2011). doi:10.1029/2010GL045225
[3] N. I. Gershenzon and G. Bambakidis, Transition from static to dynamic macroscopic friction in the framework of the Frenkel-Kontorova model, Tribology International 61, 11-18 (2013). http://dx.doi.org/10.1016/j.triboint.2012.11.025
[4] N. I. Gershenzon and G. Bambakidis, Model of deep nonvolcanic tremor part I: Ambient and triggered tremor, Bull. Seismol. Soc. Am. 104, (4) (2014). doi: 10.1785/0120130234.
[5] N. I. Gershenzon and G. Bambakidis, Model of Deep Nonvolcanic Tremor Part II: Episodic Tremor and Slip Bull. Seismol. Soc. Am., 105, (2a) (2015), doi: 10.1785/0120140225

[6] N. I. Gershenzon, G. Bambakidis and T. Skinner, How does dissipation affect the transition from static to dynamic macroscopic friction? Lubricants 3, (2) 311-331 (2015). doi:10.3390/lubricants20x000x
[7] N. I. Gershenzon, G. Bambakidis and T. Skinner, Sine-Gordon modulation solutions: Application to macroscopic non-lubricant friction, Physica D 333, 285-292 (2016).
[8] N. I. Gershenzon, Instability of frictional sliding, Bull. Seismol. Sos. Am 109 (3), 1164-1170 (2019) https://doi.org/10.1785/0120180273
[9] C. H. Scholz, The mechanics of earthquakes and faulting, University Press, Cambridge, (2002).
[10] I. S. Sacks, A. T. Linde, J. A. Snoke, and S. Suyehiro, A slow earthquake sequence following the Izu-Oshima earthquake of 1978, in Earthquake Prediction, edited by D.W. Simpson and P. G. Richards, American Geophysical Union, 617-628,Washington, D.C. (1981).
[11] T. Linde and P. G. Silver, Elevation changes and the great 1960 Chilean earthquakes: support for aseismic slip, Geophys. Res. Lett. 16, 1305-1308 (1989).
[12] H. Kanamori, A slow seismic event recorded in Pasadena. Geophys. Res. Lett. 16, (12) 1411-1414 (1989).
[13] H. Hirose, K. Hirahara, F. Kimata, N. Fujii, and S. Miyazaki, A slow thrust slip event following the two 1996 Hyuganada earthquakes beneath the Bungo Channel, Southwest Japan, Geophys. Res. Lett. 26, (21) 3237-3240 (1999).
[14] H. Dragert, K. Wang, and T. S. James, A silent slip event on the deeper Cascadia subduction interface, Science 292, (5521) 1525-1528 (2001).
[15] A. R. Lowry, K. M. Larson, V. Kostoglodov, and R. Bilham, Transient fault slip in Guerrero, SouthernMexico. Geophys. Res. Lett, 28, (19) 3753-3756 (2001). https://doi.org/10.1029/2001GL013238
[16] G. Rogers and H. Dragert, Episodic tremor and slip on the Cascadia subduction zone: The chatter of silent slip. Science 300, (5627) 1942-1943 (2003).
[17] A. Douglas, J. Beavan, L. Wallace, and J. Townend, Slow slip on the northern Hikurangi subduction interface, New Zealand, Geophys.l Res. Lett. 32, L16305 (2005). https://doi.org/10.1029/2005GL023607
[18] M. Vallee, J.-M. Nocquet, J. Battaglia, Y. Font, M. Segovia, M. Regnier, P. Mothes, P. Jarrin, D. Cisneros, S. Vaca et al., Intense interface seismicity triggered by a shallow slow slip

event in the central Ecuador subduction zone. J. Geophys. Res.: Solid Earth, 118, 2965-2981 (2013). https://doi.org/10.1002/jgrb. 50216
[19] S. Ruiz, M. Metois, A. Fuenzalida, J. Ruiz, F. Leyton, R. Grandin, C. Vigny, R. Madariaga, J. Campos, Intense foreshocks and a slow slip event preceded the 2014 Iquique Mw 8.1 earthquake. Science 345, (6201) 1165-1169 (2014). https://doi.org/10.1126/science. 1256074
[20] M. S. Boettcher and T. H. Jordan, Earthquake scaling relations for mid-ocean ridge transform faults, J. Geophys. Res., 109, B12302 (2004) doi:10.1029/2004JB003110
[21] K. Obara and A. Kato, Connecting slow earthquakes to huge earthquakes, Science 353, (6296) 253-257 (2016). https://doi.org/10.1126/science.aaf1512
[22] S. Ide, G. C. Beroza, D. R. Shelly, and T. Uchide, A scaling law for slow earthquakes, Nature 447, (7140) 76-79 (2007). https://doi.org/10.1038/nature05780
[23] H. Houston and J. E Vidale, Relationships in a slow slip, Nature 447, (7140) 49-50 (2007)
[24] H. Gao, D. A. Schmidt, and R. J. Weldon II, Scaling Relationships of Source Parameters for Slow Slip Events, Bull. Seismol. Soc. Am. 102, (1) 352-360 (2012). doi: 10.1785/0120110096
[25] Y. Ben-Zion, Collective behavior of earthquakes and faults: Continuum-discrete transitions, progressive evolutionary changes, and different dynamic regimes, Rev. Geophys., 46, RG4006, (2008), doi:10.1029/2008RG00026.
[26] B. V. Kostrov, Self-similar problems of propagation of shear cracks, J. Appl. Math. Mech. 28, 1077-1087 (1964).
[27] B. V. Kostrov, Unsteady propagation of longitudinal shear cracks, J. Appl. Math. Mech. 30, 1241-1248 (1966).
[28] B. V. Kostrov, Crack propagation at variable velocity, J. Appl. Math. Mech. 38, 511-519 (1974).
[29] R. Burridge and J. Willis, The self-similar problem of the expanding elliptical crack in an anisotropic solid, Proc. Camb. Phil. Soc. 66, 443-468 (1969).
[30] R. Burridge and C. Levy, Self-similar circular shear cracks lacking cohesion, Bull. Seismol. Soc. Am. 64, 1789-1808 (1974).
[31] R. Madariaga, Dynamics of an expanding circular fault, Bull. Seismol. Soc. Am. 66, 639666 (1976).
[32] E. Fukuyama and R. Madariaga, Rupture Dynamics of a Planar Fault in a 3D Elastic Medium: Rate- and Slip-Weakening Friction Bull. Seismol. Soc. Am. 88, (1), 1-17, (1998)

[33] D. J. Andrews, Rupture velocity of plane strain shear cracks, J. Geophys. Res. 81, 56795687 (1976).
[34] D. J. Andrews, Test of two methods for faulting in finite-difference calculations, Bull. Seismol. Soc. Am. 89, (4) 931-937 (1999).
[35] A. Bizzarri and M. Cocco, 3D dynamic simulations of spontaneous rupture propagation governed by different constitutive laws with rake rotation allowed, Ann. Geophys. 48, (2) 279299 (2005).
[36] S. Ma, S. Custo’dio, R. J. Archuleta, and P. Liu, Dynamic modeling of the 2004 Mw 6.0 Parkfield, California, earthquake, J. Geophys. Res. 113, B02301 (2008)
[37] J., Schmedes, R. J. Archuleta, and D. Lavalle’e, Correlation of earthquake source parameters inferred from dynamic rupture simulations, J. Geophys. Res. 115, B03304 (2010). doi:10.1029/2009JB006689.
[38] R. Burridge and L. Knopoff, Model and theoretical seismicity, Bull. Seismol. Soc. Am. 57, 341-371 (1967).
[39] J. Carlson and J. S. Langer, Mechanical model of an earthquake fault, Phys. Rev. A 40, 6470-6484 (1989).
[40] S. Brown, C. H. Scholz, J. B. Rundle, A simplified spring-block model of earthquakes GRL. 18(2), 215-218 (1991); doi:10.1029/2007JB005216.
[41] Huisman, B.A.H.; Fasolino, A. Transition to strictly solitary motion in the BurridgeKnopoff model of multicontact friction. Phys. Rev. E 2005, 72, doi:10.1103/PhysRevE.72.016107.
[42] O. M. Braun, I. Barel, and M. Urbakh, Dynamics of Transition from Static to Kinetic Friction, Phys. Rev. Lett. 103, 194301 (2009).
[43] J. Trømborg, J. Scheibert, D. S. Amundsen, K. Thøgersen, and A. Malthe-Sørenssen, Transition from Static to Kinetic Friction: Insights from a 2D Model, Phys. Rev. Lett. 107, 074301 (2011).
[44] J. E. Morales, G. James, A. Tonnelier, Solitary waves in the excitable Burridge-Knopoff model, Wave Motion 76, 103-121 (2018)
[45] S. T. Tse and J. R. Rice, Crustal earthquake instability in relation to the depth variation of frictional slip properties, J. Geophys. Res. 91, (B9) 9452-9472 (1986).

[46] J. H. Dieterich, Modeling of rock friction, 1. Experimental results and constitutive equations, J. Geophys. Res. 84, (B5) 2161-2168 (1979).
[47] A. L. Ruina, Slip instability and state variable friction laws, J. Geophys. Res. 88, 10359-10370 (1983).
[48] J. R. Rice, Spatio-temporal Complexity of Slip on a Fault, J. Geophys. Res. 98, (B6) 98859907 (1993).
[49] N. Lapusta, J. R. Rice, Y. Ben-Zion, and G. Zheng, Elastodynamic analysis for slow tectonic loading with spontaneous rupture episodes on faults with rate- and state-dependent friction, J. Geophys. Res. 105, 23765-23789 (2000).
[50] N. Lapusta and J. R. Rice, Nucleation and early seismic propagation of small and large events in a crustal earthquake model, J. Geophys. Res. 108, (B4) 2205 (2003). doi:10.1029/2001JB000793
[51] J.-P. Ampuero and A. M. Rubin, Earthquake nucleation on rate and state faults - Aging and slip laws, J. Geophys. Res. 113, B01302 (2008). doi:10.1029/2007JB005082
[52] Y. Liu and J. R. Rice, Aseismic slip transients emerge spontaneously in three-dimensional rate and state modeling of subduction earthquake sequences, J. Geophys. Res. 110, B08307 (2005). doi:10.1029/2004JB003424.
[53] Y. Liu and J. R. Rice, Spontaneous and triggered aseismic deformation transients in a subduction fault model, J. Geophys. Res. 112, B09404 (2007). doi:10.1029/2007JB004930
[54] B. Shibazaki and T. Shimamoto, Modelling of short-interval silent slip events in deeper subduction interfaces considering the frictional properties at the unstable-stable transition regime, Geophys. J. Int. 171, 191-205 (2007). doi:10.1111/j.1365-246X.2007.03434.x.
[55] A. M. Rubin, Episodic slow slip events and rate-and-state friction, J. Geophys. Res. 113, B11414 (2008). doi:10.1029/2008JB005642.
[56] A. M. Rubin and J.-P. Ampuero, Self-similar slip pulses during rate-and-state earthquake nucleation, J. Geophys. Res. 114, B11305 (2009). doi:10.1029/2009JB006529.
[57] P. Segall, A. M. Rubin, A. M. Bradley, and J. R. Rice, Dilatant strengthening as a mechanism for slow slip events, J. Geophys. Res. 115, B12305 (2010). doi:10.1029/2010JB007449.

[58] P. Romanet, H. S. Bhat, R. Jolivet, and R. Madariaga, Fast and slow slip events emerge due to fault geometrical complexity, Geophys Res. Lett. 45, 4809-4819 (2018). https://doi.org/10.1029/2018GL077579
[59] C. W. A. Harbord, S. B. Nielsen, N. De Paola, and R. E. Holdsworth, Earthquake nucleation on rough faults, Geology, 45(10), 931-934 (2017), doi:10.1130/G39181.1.
[60] G. G. Adams, Self-excited oscillations of two elastic half-spaces sliding with a constant coefficient of friction, J. Appl. Mech. 62 (4), 867-872, (1995). doi: 10.1115/1.2896013.
[61] H. Perfettini, and J. P. Ampuero, Dynamics of a velocity strengthening fault region: Implications for slow earthquakes and postseismic slip, J. Geophys. Res., 113(B9) (2008) doi:10.1029/2007jb005398.
[62] S. M. Rubinstein, G. Cohen, and J. Fineberg, Dynamics of precursors to frictional sliding, Phys. Rev. Lett. 98, 226103 (2007).
[63] O. Ben-David, G. Cohen, and J. Fineberg, The dynamics of the onset of frictional slip, Science 330, 211-214 (2010a). http://dx.doi.org/10.1126/science.1194777.
[64] O. Ben-David, S. M. Rubinstein, and J. Fineberg, Slip-stick and the evolution of frictional strength, Nature 463, 76-79 (2010b).
[65] D. J. Andrews, A stochastic fault model: 2. Time-dependent case. J. Geophys. Res.: Solid Earth 86, (B11) 10821-10834 (1981).
[66] J. H. Dieterich and D. E. Smith, Nonplanar faults: Mechanics of slip and off-fault damage. In Mechanics, structure and evolution of fault zones (1799-1815). Birkhäuser Basel (2009).
[67] D. E. Smith and T. H. Heaton, Models of stochastic, spatially varying stress in the crust compatible with focal-mechanism data, and how stress inversions can be biased toward the stress rate. Bull. Seismol. Soc. Am. 101, (3) 1396-1421 (2011).
[68] Z. Q. Shi and S. M. Day, Rupture dynamics and ground motion from 3-D rough-fault simulations, J. Geophys. Res. 118, (3) 1122-1141 (2013). doi:10.1002/jgrb. 50094.
[69] J. L. Jiang and N. Lapusta, Deeper penetration of large earthquakes on seismically quiescent faults, Science 352, (6291) 1293-1297 (2016) doi:10.1126/science.aaf1496.
[70] Y. Y. Lin and N. Lapusta, Microseismicity simulated on asperity-like fault patches: on scaling of seismic moment with duration and seismological estimates of stress drops, Geophys. Res. Lett. 45, (16) 8145-8155 (2018). doi:10.1029/2018gl078650.

[71] F. Yuan and V. Prakash, Slip weakening in rocks and analog materials at co-seismic slip rates, J. Mech. Phys. Solids 56, 542-560 (2008)
[72] N. M. Beeler, T. E. Tullis, and D. L. Goldsby, Constitutive relationships and physical basis of fault strength due to flash heating, J. Geophys. Res. 113, B01401 (2008). doi:10.1029/2007JB004988
[73] D. L. Goldsby and T. E. Tullis, Flash Heating Leads to Low Frictional Strength of Crustal Rocks at Earthquake Slip Rates, Science 34, (6053) 16-218 (2011) DOI: 10.1126/science. 1207902
[74] T. H. Heaton, Evidence for and implications of self-healing pulses of slip in earthquake rupture, Phys. Earth Planet. Interiors 64, 1-20 (1990).
[75] A. Cochard and R. Madariaga, Dynamic faulting under rate-dependent friction, Pure Appl. Geophys 142, 419-445 (1994).
[76] A. Cochard and R. Madariaga, Complexity of seismicity due to highly rate dependent friction, J. Geophys. Res. 101, 25321-25336 (1996).
[77] G. Perrin, J. R. Rice, and G. Zheng, Self-heating slip pulse on a frictional surface, J. Mech. Phys. Solids 43, 1461-1495 (1995).
[78] N. M. Beeler and T. E. Tullis, Self-healing slip pulse in dynamic rupture models due to velocity-dependent strength, Bull. Seism. Soc. Am. 86, 1130-1148 (1996).
[79] G. Zheng and J. R. Rice, Conditions under which velocity-weakening friction allows a selfhealing versus a cracklike mode of rupture, Bull. Seismol. Soc. Am. 88, (6) 1466-1483 (1998).
[80] D. J. Andrews and Y. Ben-Zion, Wrinkle-like slip pulse on a fault between different materials, J. Geophys. Res. 102, 553 -571 (1997) doi:10.1029/96JB02856.
[81] G. G. Adams, Steady sliding of two elastic half-spaces with friction reduction due to interface stick-slip, J. Appl. Mech. 65, 470-475 (1998). doi:10.1115/1.2789077.
[82] E. G. Daub, M. L. Manning, and J. M. Carlson, Pulse-like, crack-like, and supershear earthquake ruptures with shear strain localization, J. Geophys. Res. 115, B05311 (2010). doi:10.1029/2009JB006388
[83] A. Bizzarri, Pulse-like dynamic earthquake rupture propagation under rate-, state- and temperature-dependent friction, Geophys. Res. Lett. 37, L18307 (2010). doi:10.1029/ 2010GL044541.

[84] E.M. Dunham, D. Belanger, L. Cong, and J. E. Kozdon, Earthquake ruptures with strongly rate-weakening friction and off-fault plasticity, Part 1: Planar faults, Bull. Seismol. Soc. Am. 101, 2296-2307 (2011).
[85] Y. Huang and J.-P. Ampuero, Pulse-like ruptures induced by low-velocity fault zones, J. Geophys. Res. 116, B12307 (2011). doi:10.1029/2011JB008684.
[86] V. Lyakhovsky and Y. Ben-Zion, A continuum damage breakage faulting model accounting for solid-granular transitions, Pure Appl. Geophys., 171, 3099-3123 (2014). doi:10.1007/s00024-014-0845-4.
[87] T. A. Kontorova and Y. I. Frenkel, On the theory of plastic deformation and twinning, Zh. Eksp. Teor. Fiz. 8, 89-95 (1938).
[88] N. I. Gershenzon, Interaction of a group of dislocations within the framework of the continuum Frenkel-Kontorova model, Phys. Rev. B 50, 13308-13314 (1994).
[89] J. P. Hirth and J. Lothe, Theory of Dislocations, 2nd ed.; Krieger Publishing Company: Malabar, FL, USA, (1982).
[90] F. P. Bowden and D. Tabor, Friction and Lubrication; Oxford University Press: Oxford, UK, (1954).
[91] E. Rabinowicz, Friction and Wear of Materials, Wiley: New York, NY, USA, (1965).
[92] X-Q. Chen, H. Niu, D. Li, and Y. Li, Modeling hardness of polycrystalline materials and bulk metallic glasses, Intermetallics 19, (9), 1275-1281 (2011).
[93] O. M. Braun and Yu. S. Kivshar, The Frenkel-Kontorova Model: Concepts, Methods, and Applications, Springer-Verlag, Berlin, (2004).
[94] V.N. Nikolaevskiy, Geomechanics and Fluidodynamics: with Applications to Reservoir Engineering. Kluwer Acad Publish, Dordrecht Boston London (1996).
[95] V. G. Bykov, A model of unsteady-state slip motion on a fault in a rock sample. Izvestiya, Physics of the Solid Earth 37, (6), 484-488 (2001).
[96] V. G. Bykov, Prediction and observation of strain waves in the Earth, Geodynamics & Tectonophysics 9, (3) 721-754 (2018) doi:10.5800/GT-2018-9-3-0369.
[97] G. B. Whitham, Linear and nonlinear waves, Wiley, New York, (1974).
[98] J. R. Rice, Heating and weakening of faults during earthquake slip, J. Geophys. Res. 111, B05311 (2006) doi:10.1029/2005JB004006

[99] E. E. Brodsky, J. J. Gilchrist, A. Sagy, and C. Collettini, Faults smooth gradually as a function of slip, Earth and Planetary Science Letters 302,(1-2) 185-193 (2011). https://doi.org/10.1016/j.epsl.2010.12.010
[100] K. Obara, Nonvolcanic deep tremor associated with subduction in southwest Japan, Science 296, 1679-1681 (2002).
[101] K. Obara, Inhomogeneous distribution of deep slow earthquake activity along the strike of the subducting Philippine Sea plate, Gondwana Res. 16, 512-526 (2009). doi: 10.1016/j.gr.2009.04.011.
[102] V. Kostoglodov, S. K. Singh, J. A. Santiago, S. I. Franco, K.M. Larson, A. R. Lowry, and R. Bilham, A large silent earthquake in the Guerrero seismic gap,Mexico, Geophys. Res. Lett. 30, (15) 1807 (2003). doi: 10.1029/2003GL017219.
[103] H. Kao, S. Shan, H. Dragert, G. Rogers, J. F. Cassidy, and K. Ramachandran, A wide depth distribution of seismic tremors along the northern Cascadia margin, Nature 436, 841-844 (2005).
[104] Z. Peng and J. Gomberg, An integrative perspective of coupled seismic and aseismic slow slip phenomena, Nature Geosci. 3, 599-607 (2010). doi: 10.1038/ngeo940
[105] A. Ghosh, J. E. Vidale, J. R. Sweet, K. C. Creager, A. G. Wech, and H. Houston, Tremor bands sweep Cascadia, Geophys. Res. Lett., 37, L08301, (2010). doi:10.1029/2009GL042301
[106] E. E. Brodsky and H. Kanamori, Elastohydrodynamic lubrication of faults J. Geoph. Res. 106, B8 16357-16374 (2001)
[107] Rice, J. R. (2006), Heating and weakening of faults during earthquake slip, J. Geophys. Res., 111, B05311, doi:10.1029/2005JB004006.
[108] H. Noda, E. M. Dunham, and J. R. Rice, Earthquake ruptures with thermal weakening and the operation of major faults at low overall stress levels, J. Geophys. Res. 114, B07302 (2009) doi:10.1029/2008JB006143.
[109] M. Bouchon, M. Vallée, Observation of Long Supershear Rupture During the Magnitude 8.1 Kunlunshan Earthquake, Science 301, 824 (2003).
[110] K. Xia, A. J. Rosakis, H. Kanamori, Laboratory Earthquakes: The Sub-Rayleigh-toSupershear Rupture Transition, Science 303, (5665) 1859-1861 (2004).
DOI: 10.1126/science. 1094022

[111] S. Das and K. Aki, A numerical study of two-dimensional spontaneous rupture propagation Geophys. J. R. astr. SOC. 50, 643-668 (1977)
[112] E. M. Dunham, P. Favreau, and J. M. Carlson, A Supershear Transition Mechanism for Cracks, Science 299, (5612) 1557-1559 (2003). DOI: 10.1126/science. 1080650
[113] R. A. Harris and S. M. Day, Effects of a Low-Velocity Zone on a Dynamic Rupture Bull. Seismol. Soc. Am. 87, (5) 1267-1280 (1997)
[114] A. Bizzarri, Rupture speed and slip velocity: What can we learn from simulated earthquakes?, Earth and Planet. Sci. Lett. 317-318, 196-203 (2012).
[115] Y. Huang, W. L. Ellsworth, G. C. Beroza, Stress drops of induced and tectonic earthquakes in the central United States are indistinguishable, Sci. Adv. 3, (8) e1700772 (2017)
[116] H. Houston, B. G. Delbridge, A. G.Wech, and K. C. Creager, Rapid tremor reversals in Cascadia generated by a weakened plate interface, Nature Geosci. 4, 404-409 (2011). doi: 10.1038/NGEO1157

Tables

Table 1. Values of the parameters used in simulations

Figures

Figure 1. The steady state θ=θss\theta=\theta_{s s} is used as an example of the velocity-dependent friction coefficient, plotted as a function of (dimensionless) slip velocity, W:μ(W,θ=Dc/W)≡μss(W)W: \mu\left(W, \theta=D_{c} / W\right) \equiv \mu_{s s}(W) as given in Eq. (10) for W≤Wflash W \leq W_{\text {flash }} and μss(W)≡μ(W)\mu_{s s}(W) \equiv \mu(W) from Eq. (11) for W>Wflash W>W_{\text {flash }}, using parameters listed in Table 1 which give a velocity-weakening (a<b)(a<b) example for rate-state friction. The value μ0=μSS(W0)\mu_{0}=\mu_{S S}\left(W_{0}\right), where W0≈10−12W_{0} \approx 10^{-12} in dimensionless units has been chosen as a representative value for the average slip velocity along a fault, corresponding to ∼3 cm/yr\sim 3 \mathrm{~cm} / \mathrm{yr}.

Figure 2. Spatial distribution of initial shear stress for the three configurations given in Eqs. (12a-c), applied, without loss of generality, in the positive xx-direction. In each panel, δΣ=\delta \Sigma= Σmax⁡−Σmin⁡\Sigma_{\max }-\Sigma_{\min }, illustrated in particular for case (a), and Σs0\Sigma_{s 0} is the average shear stress value, illustrated in particular for case (b). The first two configurations mimic the spatial distribution of shear stress around an obstacle. The third configuration was chosen to represent effects due to a localized reduction of effective normal stress.

Figure 3. The spatial distribution of shear stress, ΣS\Sigma_{S}, and slip velocity, WW, are shown at three different times, as calculated from our model, Eqs. (5-12). Solutions were obtained using the initial stress distribution of Eq. 12a12 a and values listed in Table 1. The example values Σ50=0.7\Sigma_{50}=0.7 and δΣ=0.1\delta \Sigma=0.1 give results that are typical for regular EQs. The average rupture velocity, VV, is simply the location of the rupture front (as measured from the initiation point x=0x=0 ) divided by the time travelled. The region of non-zero slip velocity propagates to the left at V−=3.34 km/s=0.61csV^{-}=3.34 \mathrm{~km} / \mathrm{s}=0.61 c_{s}, and another region of non-zero slip velocity propagates to the right at V+=2.21 km/s=0.40csV^{+}=2.21 \mathrm{~km} / \mathrm{s}=0.40 c_{s}. The same holds true for regions of reduced and increased stress.

Figure 4. As in Fig. 3, showing further detail for the spatial-temporal distribution of shear stress, ΣS\Sigma_{S}, and slip velocity, WW, as well as stress drop, δΣ\delta \Sigma, and slip, SS for a regular EQ. The rupture front is clearly evident in each panel as the discontinuity between unchanged initial values, located at positions that are beyond the rupture edge, and modified values within the rupture boundary. The front, which clearly must be the same for any of the variables shown, travels in this case to the left at velocity V−V^{-}V−and to the right at V+V^{+}, as obtained in Fig. 3.

Figure 5. As in Fig. 4 (initial stress distribution of Eq. 12a), except for model conditions ΣS0=0.18\Sigma_{S 0}=0.18 and E=5⋅10−4\mathscr{E}=5 \cdot 10^{-4} that give a SSE. The rupture propagation velocities are V−≈V+≈0.88 km/V^{-} \approx V^{+} \approx 0.88 \mathrm{~km} / day.

Figure 6. The dynamics of the variables that characterize a rupture depend critically on the initiating shear to normal stress ratio, as illustrated in Figs. 3-5. In (a), the maximum (dimensionless) slip velocity, Wmax⁡W_{\max }, occurring within a rupture is plotted as a function of the peak shear stress, Σmax⁡=ΣS0+δΣ/2\Sigma_{\max }=\Sigma_{S 0}+\delta \Sigma / 2, relative to the normal stress ΣN\Sigma_{N}, for different values of δΣ\delta \Sigma. Solutions of the model Eqs. (5-12) were obtained using the initial stress distribution of Eq. 12a together with values listed in Table 1. Similar results are obtained for the distributions given in Eqs. 12 (b) and ©. The value for Wmax⁡W_{\max } can be found from the value at the rupture initiation point x=0x=0 soon after the rupture occurs and at any subsequent times (see text). Abrupt jumps from the dashed to the solid curves result from the change in rate-state friction, Eq. (10), to flash heating, Eq. (11), when W≥Wflash =1.1×10−3(0.1 m/sW \geq W_{\text {flash }}=1.1 \times 10^{-3}(0.1 \mathrm{~m} / \mathrm{s} in physical units). Since regular EQs occur only if Wmax⁡W_{\max } is larger than 1.1×10−21.1 \times 10^{-2} ( ≥1 m/s\geq 1 \mathrm{~m} / \mathrm{s} in physical units), we find that the reduction of friction due to flash heating is crucial for initiating EQ events. Otherwise, Wmax⁡W_{\max } never attains the critical value, as shown in the dotted curves in (a), which were obtained using only the ratestate friction law. Rate-state friction using experimental values for aa and bb (see Table 1), can give Wmax⁡≥10−2W_{\max } \geq 10^{-2} only for unrealistically large values of the shear to normal stress (greater than ∼0.9\sim 0.9 ). The rupture propagation velocity V−V^{-}V−is also plotted in (b), where the curves in the bottom right rectangle show the results of not including flash heating effects.

Figure 7. As in Fig. 3, except the initial stress distribution is given by Eq. (12b) (see also Fig. 2b). The example values Σ50=0.7\Sigma_{50}=0.7 and δΣ=0.17\delta \Sigma=0.17 give results that are typical for regular EQs and produce a crack-like rupture as defined in section III B and seen more clearly in Fig. 8.

Figure 8. As in Fig. 7 (initial stress distribution of Eq. 12b), showing further detail for the spatial-temporal distribution of shear stress, ΣS\Sigma_{S}, and slip velocity, WW, as well as stress drop, δΣ\delta \Sigma, and slip, SS, for a regular, crack-like rupture.

Figure 9. Similar to Fig. 7 (initial stress distribution of Eq. 12b), but using values ΣS0=0.7\Sigma_{S 0}=0.7 and δΣ=0.017\delta \Sigma=0.017 that result in a SSE. Note the reduction in δΣ\delta \Sigma compared to the value used in Fig. 7. In this case, the rupture expands in both directions by self-healing (Heaton-like) pulses.

Figure 10. As in Fig. 9, showing further detail for the spatial-temporal distribution of shear stress, ΣS\Sigma_{S}, and slip velocity, WW, as well as stress drop, δΣ\delta \Sigma, and slip, SS. The slip velocity clearly shows the pulse-like nature of the rupture.

Figure 11. Values of the average shear stress, Σ50\Sigma_{50}, and δΣ\delta \Sigma (maximum minus minimum stress), normalized to the normal stress, ΣN\Sigma_{N}, that determine whether a rupture event is crack-like or pulse-like. Simulations were performed over wide ranges of δΣ50\delta \Sigma_{50} and δΣ\delta \Sigma for the initial stress distributions given by Eqs. 12. For the localized initial stress distributions considered here, if Wmax⁡≥Wflash W_{\max } \geq W_{\text {flash }}, the rupture is crack-like, while the rupture is pulse-like for Wmax⁡<Wflash W_{\max }<W_{\text {flash }}. For a non-localized step-like stress distribution, ruptures are crack-like for all values of Σ50\Sigma_{50} and δΣ\delta \Sigma. See Eq. 12 and Fig. 2 for the stress profiles used in the simulations.

Figure 12. Spatial-temporal distribution of shear stress, ΣS\Sigma_{S}, slip velocity, WW, stress drop, δΣ\delta \Sigma, and slip, SS, for the initial stress distribution of Eq. 12c (see also Fig. 2c). The example values Σx0=0.7\Sigma_{x 0}=0.7 and δΣ=0.05\delta \Sigma=0.05 give results that are typical of regular EQs, which are characterized by crack-like ruptures. However, the position of the maximum value for SS moves to the right in time, in contrast to results for the initial distributions of Eqs. 12a12 a and 12b12 b in which the maximum remains centered at the initiation point x=0x=0.

Figure 13. As in Fig. 12, but for conditions Σ30=0.2\Sigma_{30}=0.2 and E=5⋅10−4\mathscr{E}=5 \cdot 10^{-4} that give a SSE and a pulse-like rupture. For comparison, the values of initial stress, Eq. 12b, used in Fig. 10, also produce a SSE and generate a pulse-like rupture. However, in the present case, the rupture now expands in only one direction, as shown by the spatial-temporal behavior of the slip and slip velocity. There are still two stress pulses, as for the initial condition of Eq. 12b, but the stress pulse at x<0x<0 is not accompanied by a slip pulse.