Slow-Slip Phenomena Represented by the One-Dimensional Burridge-Knopoff Model of Earthquakes

Slow-slip phenomena, including afterslips and silent earthquakes, are studied using a one-dimensional Burridge--Knopoff model that obeys the rate-and-state dependent friction law. By varying only a few model parameters, this simple model allows reproducing a variety of seismic slips within a single framework, including main shocks, precursory nucleation processes, afterslips, and silent earthquakes.

An earthquake is a stick-slip dynamic instability of a preexisting fault driven by the motion of a tectonic plate. 1) Numerical simulations of earthquakes based on a simplified statistical model, the so-called Burridge-Knopoff (BK) model, 2) are popular in statistical physics, and they provide considerable information about the statistical properties of earthquakes. 3,4) Although the BK model has been successful in describing earthquakes, almost all studies so far have been limited to the high-speed rupture of earthquakes or to main shocks.
Meanwhile, recent development in modern GPS technology and in high-density GPS and seismograph networks has revealed a rich variety of slow-slip phenomena, including afterslips, 5,6) silent earthquakes, [7][8][9][10][11][12] deep tremors, 13) etc., where the fault sliding velocity is several orders of magnitudes slower than that of the standard high-speed rupture. Thus, the concept of seismicity has been broadened dramatically. 14) Then, to gain a complete understanding of earthquake phenomena, one needs to incorporate these slow-slip phenomena.
It is a challenge to understand such a wide variety of seismicity from a general physical viewpoint, including slow slips. Therefore, questions such as what are the characteristics of slow-slip phenomena, how it differs from the standard high-speed rupture of a main shock, what conditions cause them to occur, etc., need to be answered. In the present letter, we wish to address this issue from the statistical-physics viewpoint by employing the 1D BK model obeying the rateand-state dependent friction (RSF) law. [15][16][17] We successfully reproduce a variety of seismic phenomena, including highspeed rupture of main shocks, its precursory nucleation processes, afterslips, and silent earthquakes, by varying only a few fundamental parameters of the model. In particular, regarding the occurrence of slow-slip phenomena, the relative magnitude of the frictional parameters a and b characterizing the RSF law turns out to be crucial. The 1D BK model obeying the RSF law considered here consists of a 1D array of N identical blocks with mass m, which are mutually connected by two neighboring blocks via elastic springs with spring stiffness k c , and these are connected to a moving plate via springs with spring stiffness k p . All blocks are subject to a friction force Φ. The dimensionless * E-mail:kawamura@ess.sci.osaka-u.ac.jp equation of motion for the i-th block can be written as 3,4) where ℓ ≡ (k c /k p ) 1/2 . The dimensionless displacement u i is normalized by the characteristic slip distance L associated with the RSF law, the time t by ω −1 = m/k p , the block velocity v i and the pulling speed of the plate ν by Lω, and the dimensionless friction force φ by k p L.
The RSF force φ i assumed here reads as [18][19][20] where θ i is the dimensionless state variable describing the "state" of the interface, and the normalized friction parameters a, b, and c represent velocity-strengthening, velocityweakening, and constant parts of friction, respectively. The original friction parameters A, B, and C are related to the normalized parameters by A = (k p L/N)a, B = (k p L/N)b, and C = (k p L/N)c, where N is the normal load. For simplicity, we inhibit the motion in the direction opposite to the motion of the plate. The state variable θ i is assumed to obey the aging law, 16) Note that, in contrast to the standard RSF law with the a term being proportional to a pure logarithmic form log v, which yields a pathological limit of a negatively divergent friction for v → 0, we phenomenologically introduce a modified form using the crossover velocity v * . [18][19][20] The modified form, where the a term reduces to a purely logarithmic form when v >> v * but becomes proportional to the block velocity v when v << v * , is able to describe a complete halt, unlike the standard form.
The model is characterized by only a few basic dimensionless parameters, e.g., frictional parameters a, b, and v * (c is actually irrelevant 20) ), and the elastic parameter ℓ. The plate velocity ν is also unimportant so long as it remains small, aside from setting the interseismic timescale.
Estimates of typical values of the model units representing natural faults have been given. 18) The BK model possesses a built-in time scale ω −1 corresponding to the rise time of an earthquake event. This may be estimated to be ∼ 1 [s]. The model possesses two distinct and independent length scales: one associated with the fault slip and the other with the distance along the fault. The former length scale is the critical slip distance L, which was estimated to be ≃ 1 [cm], while the other is the distance the rupture propagates per unit time, v s /ω, which was estimated to be ≃ 2 − 3 [km], where v s is the s-wave velocity along the fault. The typical plate velocity is several [cm/year] and corresponds to a very small number of ν ∼ 10 −7 − 10 −8 . The spring constant k p was related to the normal stress as N k p L ≃ 10 2 − 10 3 . Then, as C is known to take a value around 2 3 , c would be of order 10 2 -10 3 , with a and b being one or two orders of magnitude smaller than c. The crossover velocity v * is hard to estimate, though it should be much smaller than unity, and we take it as a parameter.
The properties of the high-speed rupture of the model were investigated, [20][21][22] together with the properties of its precursory nucleation process. 18,19) The two regimes with mutually different seismic properties are identified [18][19][20] as the weak and strong frictional instability regimes. In the former, the frictional instability parameter b is smaller than a critical value b c = 2ℓ 2 + 1 determined by the elastic parameter ℓ, and the main shocks possess eminently "characteristic" features 4,20) accompanied by a nucleation process with a quasi-static, long-lasting initial phase. In the latter, b > b c , and main shocks possess more or less "critical" features 4,20) with a nucleation process that is unaccompanied by the initial phase. These earlier studies on the model concentrated on the parameter regime where the high-speed rupture dominates, which corresponded to a < b, i.e., the frictional weakening dominates over strengthening. In fact, in such a "main-shock regime", the parameters a and v * turned out to be irrelevant.
Generally, since the frictional parameters a and b compete in their function, either a < b or a > b might significantly affect the model dynamics. More specifically, one may expect that, when a > ∼ b, the compensation effect due to the a term might induce slow-slip phenomena in the model. In the present letter, we perform a systematic survey of the parameter range of a > ∼ b, where the velocity strengthening a term is expected to play a significant role. We find that slow slips actually come into play there. When a is comparable to or only a few times larger than b, say, b < ∼ a < ∼ 2b, the model still exhibits a main shock of high-speed rupture, but subsequently exhibits a slow and long-lasting afterslip. We call this parameter regime "region II" in distinction with "region I" of highspeed rupture where a < ∼ b. When a is considerably larger than b, say, a > ∼ 2b, the model no longer exhibits high-speed rupture, but instead exhibits a long-lasting slow slip only, i.e., a silent earthquake. We call this parameter regime "region III".
In our simulations, for the sake of computational feasibility, we concentrate on the strong frictional instability regime The parameters a and v * are varied in the range of 1 − 240 and 10 −1 − 10 −3 , respectively, while fixing ν = 10 −7 . The system size (total number of blocks) is 800, and open boundary conditions are imposed. Concerning the initial conditions, all blocks are assumed to be at rest, i.e., v i = 0 (1 ≤ i ≤ N) at t = 0, the state variable is taken to be uniform θ i = 10 8 , while the displacement of each block is assumed to take random values uniformly distributed between -0.5 and 0.5 from block to block. Events at earlier times are transient and are affected in a non-stationary manner by the initial conditions. We For events accompanying slow slips, we show in Fig.2 the time t dependence of the sliding velocity v on a log-log plot, for a typical block located in the rupture-propagating part. The data for v * = 10 −3 are also shown together with those for v * = 10 −2 . For events belonging to region II, the maximum velocity of v max ≃ 10 1 − 10 3 is reached and is characteristic of a high-speed rupture, which is followed by long-lasting slow afterslips. In contrast, v max of the events belonging to region III is of order 10 −1 − 10 −2 , orders of magnitudes smaller than that of a high-speed rupture. Hence, the events in region III might well be identified as silent earthquakes. Its duration time is longer than that of a main shock in region I by a factor of approximately 10 4 − 10 5 for v * = 10 −3 . As v gets smaller beyond v * , the block comes to a stop fairly soon. Namely, v * signals the end of the slow-slip process. Then, a longer-lasting slow slip is naturally expected for a smaller v * .
One might wonder if the distinction between regions II and III were only gradual. We show in Fig.3 the a-dependence of v max spanning both regions II and III for various v * -values. As can be seen from the figure, a clear change of behavior occurs around a ≃ 50 − 60, although the borderline value seems slightly v * -dependent. A similar behavior is observed in the a-dependence of the rupture-propagating velocity v r shown in the inset, which corresponds to the inverse slope of the diagonal lines of Fig.1. We regard the observed change of behavior, which tends to be more eminent for smaller v * , as a signature discriminating regions II and III. Presumably, while the change is a crossover for a finite v * -value, it may become a sharp "transition" in the v * → 0 limit.
Regions I and II are discriminated by the absence/presence of an afterslip following the high-speed rupture of a main shock. This change turns out to occur roughly around a ∼ b. In the insets of Figs.1(a)-(c), we show typical event sequences in the stationary state, in which the displacement of each block is plotted each time all blocks are at rest. As can be seen from Fig.1(a), all events in region I are "partial" breaking in particular parts of the entire system with nontrivial size distribution. By contrast, as can be seen from Fig.1(c), events in region III are "entire" breaking the entire system that repeat quasi-periodically. The situation in region II is a bit more complex. In the parameter region close to region III, events tend to be "entire", as shown in the inset of Fig.1(b). Meanwhile, in the parameter region close to region I, the system is often divided into several large sectors, a situation somewhat intermediate between "partial" (region I) and "entire" (region III) ruptures.
In the main shock of region I, the stress distribution after the event tends to be spatially inhomogeneous, which serves as a "barrier" for subsequent events. Meanwhile, in the events in regions II and III, slow slips tend to smoothen such spatial inhomogeneity, eliminating the potential barrier for subsequent events.
Slips with low sliding velocity are also realized in the nucleation process preceding main shocks. 18,19) In fact, they are also realized even in silent earthquakes of region III. How such slow slips associated with the nucleation process resemble or differ from the slow slips studied here, and how they connect to the latter, are interesting questions. To gain insights into the issue, we show in Fig.4 the time dependence of several quantities for an epicenter block which experienced a slow nucleation process, including the sliding velocity v, the state variable θ, and their multiple vθ for the events in regions I (a), II (b), and III (c), each corresponding to the events shown in Figs.1(a)-(c).
As can be seen from Figs.4(b) and 4(c), the relation vθ = 1 is nearly satisfied in the slow-slip regimes. In fact, this relation implies an approximate stationarity, since θ is nearly constant from Eq.(3), and v ≃ 1/θ. Thus, if v takes a small value (or θ takes a large value) at the point where the relation vθ ≃ 1 is met, long-lasting slow slips might follow, which actually seems to happen in the afterslip in region II, and in the silent earthquake in region III. Even in region I, the relation vθ = 1 is met, but only during the short period of high-speed rupture of a main shock around its maximum v. 19) The observation that the relation vθ ≃ 1 is met during the slow slips of afterslips and silent earthquakes indicates that Letter Author Name they are fundamentally different in nature from those of the precursory nucleation process where vθ >> 1, as can be seen from Fig.4. During the slow nucleation process, θ is of order 10 10 , representing a strongly stuck interface formed during the long interseismic period, and vθ is also large of order 10 7 , in spite of its low v-value. Then, as is evident from Eq.(3), θ actually decreases rapidly from its extremely large initial value. As such, the slow slip in the nucleation process is a far-from-stationary phenomenon in spite of its low sliding velocity. By contrast, at the slow slips realized in afterslips and silent earthquakes, the strongly stuck interface has been more or less released, with θ reduced to moderate values of order 10 0 ∼ 10 4 . The relation vθ ≃ 1 is then well met, leading to near-stationarity. Another point to be noticed might be that the slow slip at silent earthquakes propagates in space with a constant rupture velocity, as can be seen from Fig.1(c), whereas the slow slip in the precursory nucleation process is spatially localized into the compact "seed" area. Note that the relation vθ = 1 during the slow-slip process does not hold strictly, but only approximately. In fact, vθ is slightly smaller than unity throughout the slow-slip process, leading to a gradual increase of θ and an associated decrease of v. Finally, when v reduces below v * , the a term changes its form due to the crossover effect, and the block motion comes to a complete stop, as shown in the insets of Figs.4(b) and (c). In other words, slow slip phenomena in afterslips and silent earthquakes are near-stationary, but are still essentially nonstationary phenomena.
Finally, to examine the relevance of our present results to real seismic observations, we attempt order estimates of the slow-slip duration time T . Fig.2 indicates that the local duration time (the duration time for a given block) of a silent earthquake in region III might be T ≃ 3 × 10 5 for v * = 10 −3 , which reads as several days in real time. Such a time scale may correspond to a short-term slow slip in real seismicity. 23,24) In fact, in our simulations, T turns out to be inversely proportional to the largely unknown model parameter v * in the v * -range studied (10 −1 ≥ v * ≥ 10 −3 ). If one extrapolates this to still smaller values of v * , one gets T ≃ 1 year for v * = 10 −5 , which may correspond to a long-term slow slip. 7,23,24) In summary, we studied the seismic properties of the 1D BK model obeying the RSF law over a broad parameter space, and we found a variety of slow-slip phenomena, including afterslips and silent earthquakes, which resembled those observed in real seismicity. We clarified the occurrence condition of such slow-slip phenomena, and the essential difference from the slow seismic slips as realized in the precursory nucleation process. Interestingly, the 1D BK model obeying the RSF law, being quite a simple model, is capable of reproducing rich seismic phenomena such as a main shock, a precursory nucleation process, an aftershock, and a silent earthquake by varying only a few model parameters.