Dynamics of earthquake nucleation process represented by the Burridge-Knopoff model

Dynamics of earthquake nucleation process represented by the Burridge-Knopoff model

Yushi Ueda    Shouji Morimoto    Shingo Kakui    Takumi Yamamoto    Hikaru Kawamura kawamura@ess.sci.osaka-u.ac.jp Graduate School of Science, Osaka University, Toyonaka, Osaka 560-0043, Japan
July 24, 2019
Abstract

Dynamics of earthquake nucleation process is studied on the basis of the one-dimensional Burridge-Knopoff (BK) model obeying the rate- and state-dependent friction (RSF) law. We investigate the properties of the model at each stage of the nucleation process, including the quasi-static initial phase, the unstable acceleration phase and the high-speed rupture phase or a mainshock. Two kinds of nucleation lengths and are identified and investigated. The nucleation length and the initial phase exist only for a weak frictional instability regime, while the nucleation length and the acceleration phase exist for both weak and strong instability regimes. Both and are found to be determined by the model parameters, the frictional weakening parameter and the elastic stiffness parameter, hardly dependent on the size of an ensuing mainshock. The sliding velocity is extremely slow in the initial phase up to , of order the pulling speed of the plate, while it reaches a detectable level at a certain stage of the acceleration phase. The continuum limits of the results are discussed. The continuum limit of the BK model lies in the weak frictional instability regime so that a mature homogeneous fault under the RSF law always accompanies the quasi-static nucleation process. Duration times of each stage of the nucleation process are examined. The relation to the elastic continuum model and implications to real seismicity are discussed.

I I. Introduction

Owing to the recent development of the GPS technology, it has been recognized now that general forms of seismic activity could be of a rich variety, often including various types of slow slip events, e.g., a preslip, an afterslip, a slow earthquake, etc. From the viewpoint of earthquake forecast, a preslip, i.e., a slow-slip event occurring prior to a mainshock, would be of special significance. Such a preslip prior to a mainshock is usually associated with the nucleation process which occurs preceeding the high-speed rupture of the mainshock. Then, a wide-spread expectation is that a large earthquake might be preceded by a precursory nucleation process which occurs prior to the high-speed rupture of the mainshock. Such a precursory phenomenon preceding mainshocks, if any, would be of paramount importance in its own right as well as in its possible connection to earthquake forecast.

Although the nucleation phenomena preceding the main rupture have been more or less confirmed by laboratory rock experiments Latour (); McLasky (), its nature, or even its very existence, remains less clear for real earthquakes Ando (); Ohta (); Kato (); Bouchon (); Tape (). We note that a similar nucleation process is ubiquitously observed in various types of failure processes in material science and in engineering.

Nucleation process is supposed to be localized to a compact “seed” area with its rupture velocity orders of magnitude lower than the seismic wave velocity Dieterich92 (); Ohnaka00 (); Ohnaka03 (); Scholzbook (); Dieterich09 (). The fault spends a very long time in this nucleation process, and then at some point, exhibits a rapid acceleration process accompanied by a rapid expansion of the rupture zone, finally getting into the final high-speed rupture of a mainshock.

The earthquake nucleation process might proceed via several distinct steps or “phases”. According to Ohnaka, it starts with an initial quasi-static process, and gets into the acceleration phase when the nucleus diameter exceeds a nucleation length , where the system gets out of equilibrium, and rapidly increases its slip velocity Ohnaka00 (); Ohnaka03 (). Then, when the nucleus diameter exceeds another nucleation length , the fault eventually exhibits a high-speed rupture of a mainshock. In this picture, there appear two characteristic length scales for the nucleus, and . These two nucleation lengths divide the nucleation process into “the initial phase” in which the nucleus size is smaller the (), “the acceleration phase” in which the nucleus size exceeds but is still smaller than (), and “the high-speed rupture phase” of mainshock ().

Because of a slow character of the slip, earthquake nucleation process might also be regarded as a type of more general slow-slip phenomena, including afterslips and slow earthquakes. While the relation between these different types of slow seismic processes poses an interesting and important question, we focus in the present paper on the nucleation process prior to the high-speed rupture of mainshock.

Under such circumstances, in order to get deeper understanding of the physical process behind the seismic nucleation process and the subsequent mainshock, a theoretical or a numerical study based on an appropriate model of an earthquake fault would be important and helpful. In such a physical modeling, the friction force is a crucially important part. The friction force now standard in seismology is the so-called rate and state dependent friction (RSF) law Dieterich79 (); Ruina (); Marone (); Scholz98 (). In the pioneering study, Dieterich derived a formula describing the nucleation length based on such RSF law Dieterich92 (). The most standard form of the nucleation length based on the RSF law might be

(1)

where is the normal stress, is the rigidity, is a characteristic slip distance and is a constant, while and are the frictional parameters associated with the RSF law, each representing the velocity-strengthening and the frictional-weakening parts of the friction. Dieterich also suggested that under certain conditions the nucleation length might be given by Dieterich92 ()

(2)

in which the parameter did not appear.

The RSF law has been used in many of numerical simulations on earthquakes, mostly in the continuum model TseRice (); Stuart (); Horowitz (); Rice (); Ben-ZionRice (); KatoHirasawa (); Bizzarri06a (); Bizzarri06b (), including the earthquake nucleation process. In particular, Ampuero and Rubin studied the properties of the nucleation process for the continuum model under the RSF law, with two representative evolution laws, i.e., the aging law AmpueroRubin () and the slip law RubinAmpuero () within the quasi-static approximation neglecting the inertia effect.

Meanwhile, a further simplified discrete model has also been used in earthquake studies. Particularly popular is the so-called spring-block model or the Burridge-Knopoff (BK) model BK (), in which an earthquake fault is modeled as an assembly of blocks mutually connected via elastic springs which are subject to the friction force and are slowly driven by an external force mimicking the plate drive.

The model might be better justified in the situation where there exists a well-developed mature fault layer, presumably corresponding to the low-velocity fault zone Huang () observed in many mature faults Ueda (). The fault layer is supposed to be uniformly pulled by the more or less rigid crust contingent to it. Because of its simplicity, the BK model is particularly suited to the study of statistical properties of earthquakes, since it often enables one to generate sufficiently many events, say, hundreds of thousands of events, to reliably evaluate its statistical properties. This type of model might also be relevant to the description of other stick-slip-type phenomena such as landslides Viesca ().

In many numerical simulations of the BK model, while a simple velocity-weakening friction law in which the friction force is assumed to be a single-valued decreasing function of the velocity has often been used CL89a (); CL89b (); CLST (); Carlson91 (); Carlson91b (); CLS-review (); Schmittbuhl (); MoriKawamura05 (); MoriKawamura06 (); MoriKawamura08a (); MoriKawamura08b (); MoriKawamura08c (); Kawamura-review (), a more realistic RSF law was also employed in some of recent numerical simulations of the model. For example, Cao and Aki performed a numerical simulation by combining the one-dimensional (1D) BK model with the RSF law in which various constitutive parameters were set nonuniform over blocks CaoAki (). Ohmura and Kawamura extended an earlier calculation by Cao and Aki to study the statistical properties of the 1D BK model combined with the RSF law with uniform constitutive parameters OhmuraKawamura (); Kawamura-review (). Clancy and Corcoran also performed a simulation of the model based on a modified version of the RSF law Clancy09 ().

Of course, the space discretization in the form of blocks is a crude approximation to the original continuum crust. It introduces the short-length cut-off scale into the problem in the form of the block size, which could in principle give rise to an artificial effect not realized in the continuum. Indeed, such a criticism against the BK model was made in the past Rice (). Rice criticized that the discrete BK model with the simple velocity-weakening law was “intrinsically discrete”, lacking in a well-defined continuum limit, arguing that the spatiotemporal complexity observed in the discrete BK model was due to an inherent discreteness of the model, which should disappear in continuum Rice (). In contrast to the simple velocity-weakening law, the RSF law possesses an intrinsic length scale corresponding the characteristic slip distance . If the grid spacing was taken smaller than the nucleation length which is proportional to this characteristic slip distance , the system tended to exhibit a quasi-periodic recurrence of large events, whereas, if the grid spacing was taken larger than it, the system exhibited an apparently complex or critical behavior. Note that the extent of the discreteness may be regarded as a measure of the underlying spatial inhomogeneity Rice ().

This problem of the continuum limit of the BK model was also addressed within the velocity-weakening friction law by Myers and Langer Myers93 (), by Shaw Shaw94 (), and by Mori and Kawamura MoriKawamura08c (), where the Kelvin viscosity term was introduced to produce a small length scale allowing for a sensible continuum limit.

In fact, the problem of the small length scale of the BK model is closely related to the nucleation phenomena. According to Rice, the continuum system under the RSF law always exhibits a quasi-static nucleation process prior to a mainshock Rice (). In view of such a claim, and also of the general importance of the seismic nucleation phenomena, in the present paper we wish to clarify by means of extensive numerical simulations on the 1D BK model how the nucleation process of the discrete BK model behaves in its continuum limit, by systematically varying the extent of the discreteness of the model. By so doing, we also wish to clarify the nature of the nucleation process of the model, and discuss its implications to real seismicity.

Of course, since the real fault plane is more or less two-dimensional (2D), the assumed one-dimensinality of our present model is a high simplification. While earlier studies suggested that most of qualitative features of the mainshock statistical properties could be captured even by the 1D model, there exists a possibility that, concerning the nucleation properties, the 2D model might exhibit a behavior different from the 1D model due to the richness of the underlying geometry Kawamura-review (). In the present paper, leaving a systematic study of the nucleation process of the 2D BK model as a future task, we first study the 1D model which is far simpler than the corresponding 2D model.

In such a way, on the basis of the 1D BK model, we wish to shed light on the nucleation process of a mature fault, e.g., the nucleation dynamics, the nucleation lengths and the duration times of each phase of the nucleation process. How these quantities depend on material parameters, and are related or unrelated to the size of the ensuing mainshock ? Such an issue would be of special significance from the standpoint of utilizing earthquake nucleation phenomena for a possible earthquake forecast. For example, if the nucleation length or is correlated with the mainshock size, e.g., a larger earthquake for a larger or , one might have a chance to predict the size of the mainshock from the measurement of the nucleation lengths. If, on the other hand, the nucleation length or is not correlated with the mainshock size, the prediction of the mainshock size from the measurement of the nucleation lengths would be impossible.

By its nature, the fault sliding velocity in the nucleation process tends to be very low. Hence, for any practical detection, it would crucially be important to clarify how fast the fault sliding velocity is, and how much time is left before the ensuing mainshock. With these motivations in mind, we try to conduct a systematic numerical and analytic study of the BK model in the following part of the paper.

A preliminary account of our simulations was already given in Ueda (). In the present paper, we conduct a systematic survey of more general parameter space, give detailed analytical treatments, and make comparison between theoretical and numerical results.

The rest of the paper is organized as follows. In section II, we define our model, the 1D BK model obeying the RSF law, and present its equation of motion. Its continuum limit is also given. In section III, we report on the results of our numerical simulations on the dynamics of the nucleation process of the model. In subsection III[A], we first illustrate mains features of the nucleation process. Two distinct parameter regimes exist, i.e., the weak frictional instability regime and the strong frictional instability regime. The two kinds of nucleation lengths and are identified. In the subsequent subsections [B]-[D], we present our numerical data on the dynamics of the model at each stage of its nucleation process, i.e., [B] the initial phase of the weak frictional instability regime, [C] the acceleration phase of the weak frictional instability regime, and [D] the acceleration phase of the strong frictional instability regime. In section IV, we report on the results of our theoretical analyses of the nucleation process of the model. After explaining in subsection [A] the basic scheme of the perturbation method employed, we examine the dynamics of the model in some detail in the following subsections [B]-[D], i.e., [B] the initial phase, [C] the acceleration phase at which the epicenter-block sliding velocity is smaller than the crossover velocity , , and [D] the acceleration phase at . Analytic expressions of the nucleation length and of the condition discriminating the weak and the strong frictional instability regimes is derived in subsection [B]. In subsection [E], we perform a mechanical stability analysis to re-derive and the weak/strong instability condition, which confirms the results from the perturbation analysis. In section V, we present the results of our numerical simulations focusing on various statistical properties characterizing the nucleation process, including the nucleation lengths and , and the duration times of each phase, averaged over many events. Their continuum limits are also examined. Finally, section VI is devoted to summary and discussion. Implications to real seismicity and possible extensions of the present analysis are discussed.

Ii II. The model and its continuum limit

The 1D BK model consists of a 1D array of identical blocks of the mass , which are mutually connected with the two neighboring blocks via the elastic springs of the spring stiffness , also connected to the moving plate via the springs of the spring stiffness , and are driven with a constant rate . All blocks are subject to the friction force , which is the source of the nonlinearity in the model.

The equation of motion for the -th block can be written as

(3)

where is the time, is the displacement of the -th block, and is the friction force at the -th block. For simplicity, the motion in the direction opposite to the plate drive is inhibited by imposing an infinitely large friction for .

For the friction law, we assume the RSF law given by

(4)

where is the sliding velocity of the -th block, is the time-dependent state variable (with the dimension of the time) representing the “state” of the slip interface, is a crossover velocity underlying the RSF law, is an effective normal load, is a critical slip distance which is a measure of the sliding distance necessary for the surface to evolve to a new state, with and positive constants describing the RSF law. The first term (-term) is a constant taking a value around , which dominates the total friction in magnitude, the second term (-term) a velocity-strengthening direct term describing the part of the friction responding immediately to the velocity change, the third part (-term) an indirect frictional-weakening term dependent of the state variable . Laboratory experiments suggest that the - and -terms are smaller than the -term by one or two orders of magnitudes, yet they play an essential role in stick-slip dynamics Marone (); Scholz98 (); Scholzbook ().

Note that, in the standard RSF law, the -term is often assumed to be proportional to . Obviously, this form becomes pathological in the limit because it gives a negatively divergent friction. In other words, the pure logarithmic form of the -term cannot describe the state at rest. We cure this pathology by phenomenologically introducing a modified form given above. The modified form, where the -term becomes proportional to the block velocity at but reduces to the purely logarithmic form at , enables one to describe a complete halt. The characteristic velocity represents a crossover velocity, describing the low-velocity cutoff of the logarithmic behavior of the friction.

For the evolution law of the state variable , we use here the so-called aging (slowness) law given by

(5)

Under this evolution law, the state variable grows linearly with the time at a complete halt , reaching a very large value at the outset of the nucleation process, while it decays very rapidly during the seismic rupture.

The equation of motion can be made dimensionless by taking the length unit to be the critical slip distance , the time unit to be and the velocity unit to be ,

(6)
(7)

where the dimensionless variables are defined by , , , , , , , , , while is the dimensionless elastic stiffness parameter.

It is sometimes more convenient to rewrite the equation of motion in terms of the velocity variable instead of the displacement . By differentiating (6) with respect to and by using (7), one gets

(8)

The block displacement can be obtained up to a constant by integrating the velocity with respect to .

One sees from eqs.(8) and (7) that the constant frictional parameter no longer remains in the governing equations, meaning this parameter is essentially irrelevant to the dynamical properties of the model. In our simulations, we use either eq.(6) or (8) depending on the situation. In solving the high-speed motion, we use eq.(6), while in solving the low-speed motion as realized in the initial phase or the early stage of the acceleration phase, we use eq.(8).

The frictional parameter / tends to suppress/enhance the frictional instability. The earthquake instability is driven primarily by the velocity-weakening -term, while the velocity strengthening -term tends to mitigate the unstable slip toward the aseismic slip. Since the frictional parameters and compete in their functions, either or might affect the dynamics significantly. Earthquake properties in this regime of will be reported in a separate paper, with emphasis on the slow-slip phenomena intrinsic to this regime. Meanwhile, we find that the properties of the precursory nucleation process of the model, which occurs preceding a mainshock, do not much depend on the relative magnitude of and . Although we study in the present paper the nucleation process in the parameter range where the unstable seismic character is dominant in a mainshock, main qualitative features of the nucleation process would not change much even for .

The setting assumed in the BK model in terms of an earthquake fault embedded in the 3D continuum crust was examined in Ueda (). In particular, estimates of typical values of the model parameters for natural earthquake faults are given as [s], a few [cm], and . This yields , around while and being one or two orders of magnitude smaller than . The crossover velocity and its dimensionless counterpart is hard to estimate though it should be much smaller than unity, and we take it as a parameter in our simulations.

The continuum limit of the BK model corresponds to making the dimensionless block size , defined by , to be infinitesimal , simultaneously making the system infinitely rigid with MoriKawamura08c (). The dimensionless distance between the block and is given by

(9)

As discussed in MoriKawamura08c (), the 1D equation of motion in the continuum limit is given in the dimensionful form by

(10)

where is the displacement at the position and the time , is the friction force per unit mass, while and are the characteristic frequency and the characteristic wave-velocity (-wave velocity), respectively. Note that the term representing the plate drive is absent in the standard elasto-dynamic equation.

Iii III. Simulation results I

In this section, we present the results of our numerical simulations on the dynamical properties of the model. After surveying their main features in subsection [A], we present detailed data in the following subsections separately for each phase of the nucleation process.

Figure 1: Color plots of typical earthquake nucleation processes depicted in the block-number (position) versus the time plane, realized in the weak frictional instability regime. The color represents the block sliding velocity (white means exactly zero slip, while a low non-zero velocity is represented by black). The model parameters are , , , , and . and are two types of nucleation lengths. Events are taken from the event sequence occurring in the stationary state of the model.

iii.1 A. Weak versus strong frictional instability regimes

The first question might be whether the 1D BK model under the RSF law ever exhibits a nucleation process prior to a mainshock, and if it does, under what conditions. Remember that our constitutive law allows for a complete stick (i.e., for all ) during the interseismic period, which enables us to unambiguously define the onset of the nucleation process by the point where one of the blocks gains a nonzero velocity.

Figure 2: The color plots representing the evolutions of the rupture process when the external loading is artificially stopped, (a) at a point when the number of moving blocks (rupture-zone size) is less than , i.e., , or (b) at a point when the number of moving blocks is greater than , i.e., , in the case of the weak frictional instability. The time origin is set to the point of here. The color represents the block sliding velocity (white means exactly zero slip). The model parameters are , , , , and , corresponding to .

We illustrate in Fig.1 a typical example of seismic events realized in the 1D BK model under the RSF law, where the time evolution of the movement of each block is shown as a color plot. In the figure, the two types of nucleation lengths and are also illustrated. The nucleation process including the initial phase and the accerrelation phase is realized in addition to the high-speed rupture phase (mainshock). The model parameters are set to , , , , and . The origin of the time () is taken to be the onset of the nucleation process where an epicenter block begins to move. An example shown in Fig.1 is an event occurring in the stationary state of the seismic sequence of the model, realized after transient initial events where the memory of initial conditions are still remnant.

As can clearly be seen from the figure, a slow nucleation process with a long duration time of order is observed. Among the two types of nucleation lengths, and () in Fig.1, is the length separating stable and unstable ruptures. Namely, when the nucleus size is less than , the rupture process is stable and reversible, whereas, when exceeds , it becomes unstable and irreversible.

One illustrative way to demonstrate the expected borderline behavior across the nucleation length may be to artificially stop the external loading in the course of a simulation. Indeed, when the external loading is stopped at a point before , the rupture itself also stops there, as demonstrated in Fig.2(a), whereas, if the external loading is stopped at any point beyond , the subsequent seismic rupture is no longer stoppable and evolves until its very end, as demonstrated in Fig.2(b). (Even better criterion might be whether the block sliding velocity is increased or decreased when the loading is artificially stopped, rather than whether the block is completely stopped or not.)

Ohnaka suggested that, in addition to the nucleation length , there exists another nucleation length , which discriminates between the acceleration phase and the high-speed rupture phase Ohnaka00 (); Ohnaka03 (). In the high-speed rupture phase beyond , the rupture propagates with a nearly constant speed in both directions in the form of two separate packets of moving blocks, as can be seen in Fig.1. In the figure, the high-speed rupture of a mainshock corresponds to the linear portion of the rupture propagation line with its slope being the propagation speed of .

While there might be several ways to define the nucleation length (), we tentatively give one definition here. As can be seen from Fig.1, the number of simultaneously moving blocks tends to be large around . Hence, we tentatively define by the size of the nucleus at which the number of simultaneously moving blocks becomes maximum for a given event, which we denote : See Fig.10(a) below. At or very close to this point, the epicenter block ceases to move and the group of moving blocks are detached into two parts, each part propagating in opposite directions. Below in subsection D, we shall give another, perhaps physically more appropriate definition of , which is actually the one indicated as in Figs.1 and 3.

In fact, we find that the slow and reversible nucleaton process corresponding to the initial phase is realized only in the “weak frictional instability” regime where the frictional-weakening parameter is smaller than a critical value , while, in the “strong frictional instability” regime where the frictional-weakening parameter is greater than a critical value , the nucleation process does not accompany the slow and reversible initial phase. A typical example of seismic events in the strong frictional instability regime is given in Fig.3, where the model parameters are set to and , other parameters being common with those of Fig.1. An apparently nucleation-like process seen in Fig.3 just before the high-speed rupture propagation is not a quasi-static initial phase, but is an unstable acceleration phase. Its duration time is around which is by many orders of magnitude shorter than the duration time of the quasi-static initial phase seen in Fig.1, . The acceleration phase in the stronger frictional instability regime sometimes could be longer, say, , particularly for smaller . Yet, the dynamics is already irreversible there.

Figure 3: Color plots of typical earthquake nucleation processes depicted in the block-number (position) versus the time plane, realized in the strong frictional instability regime. The color represents the block sliding velocity (white means exactly zero slip). The model parameters are , , , , and (, , and are taken to be common with those in Fig.1). Events are taken from the event sequence occurring in the stationary state of the model.

The borderline between the weak and the strong frictional instability regimes is given by a critical value of , , which is found to have a simple expression . This expression was already reported in Ref.Ueda (), but we shall derive this expression in several ways in the following part of the paper.

We note that large events in the weak frictional instability regime always accompany precursory nucleation process irrespective of each individual event, or the choice of the initial conditions, while, those in the strong frictional instability does not accompany a precursory nucleation process in any condition. Hence, for a given set of model parameters, the presence or absence of the quasi-static nucleation process is uniquely determined, not depending on each individual event, but is determined simply by the condition of the friction parameter being either greater or smaller than the critical value , being the elastic stiffness parameter.

In order to demonstrate the spatiotemporal evolution of the nucleation process of the model, we show in Figs.4 and 5 the time evolutions of the spatial profile of (a) the block sliding velocity , (b) the state variable , and (c) the multiple of the two quantities , in a typical nucleation process of a large event realized in the stationary state in the weak frictional instability regime. Fig.4 covers the time regime from the onset of the nucleation process till the system reaches , whereas Fig.5 from the point of till an earlier stage of the high-speed rupture phase. The model parameters are set to and . From these figures, the manner how the nucleus grows and how the nucleation process transforms into the high-speed rupture of a mainshock is clearly visible.

Some characteristic points of the nucleation process corresponding to and (in Fig.4), (in Figs.4 and 5), and (in Fig.5) are indicated by blue curves. Here is a characteristic crossover velocity at which the inertia effect becomes significant, to be defined below in §IVD. Note that, beyond the point , the inertia effect plays an important role, and the quasi-static approximation in no longer valid. As such, the time range beyond is not covered by RubinAmpuero () who employed the quasi-static approximation. In the range up to , the profiles obtained here look similar to the ones given in RubinAmpuero () for the continuum model.

Figure 4: The time evolutions of the spatial profile of (a) the block sliding velocity , (b) the state variable , and (c) the multiple of the two quantities , during a typical nucleation process in the weak frictional instability regime. The model parameters are set to and . The epicenter block is located at the center, . The covered time range is from the onset of the nucleation process until the point of . The blue curves represent the points of , and : See the text for the definitions of , and .

As can be seen from Fig.4(a), the sliding velocity gets larger until . Beyond this point, first the epicenter block, and subsequently the neighboring blocks, begin to decelerate, and eventually come to stop (Fig.5(a)). The nucleus is detached into two parts, each of which propagates in the opposite directions forming a rupture front of a mainshock.

As can be seen from Fig.4(b), in an earlier period of the nucleation process, the state variable maintains its large value acquired during the halt period between mainshocks, while it rapidly decreases in the later period as the block movement accelerates, and eventually reaches a minimum value around , first at the epicenter block, and subsequently at the neighboring blocks. After this point, the -value tends to be recovered again (Fig.5(b)).

The multiple of and , , plays an important role in the healing process since it appears on the r.h.s. of the equation of motion of the state variable, eq.(7). As can be seen from Figs.4(c) and 5(c), this quantity tends to increase in the earlier period of the nucleation process, first gradually and more rapidly beyond , reaches a maximum at a point between and , then drops very sharply until it tends to stay around a value close to unity. Note that is a special point corresponding to the stationary condition for the time evolution of the state variable: see eq.(7). Such a plateau-like behavior of arises around in the epicenter region, and transmits outwards in the nucleus. Further beyond , tends to decrease again, first in the epicenter region, and subsequently in the outer region in the nucleus.

Figure 5: The time evolutions of the spatial profile of (a) the block sliding velocity , (b) the state variable , and (c) the multiple of the two quantities , during a typical nucleation process in the weak frictional instability regime. The model parameters are set to and . The epicenter block is located at the center, . The covered time range is from the point of until the earlier stage of the high-speed rupture of a mainshock, following the time range covered by Fig.4. The blue curves represent the points of and : See the text for the definitions of and .

In the following subsections, we present our simulation data in some detail in each phase of the nucleation process, i.e., (B) the initial phase of the weak frictional instability regime, (C) the acceleration phase of the weak frictional instability regime, and (D) the acceleration phase of the strong frictional instability regime, respectively.

iii.2 B. The initial phase of the weak frictional instability regime

Let us begin with the nucleation process in the weak frictional instability regime at .

Figure 6: The time evolutions of (a) the sliding velocity , (b) the state variable , and (c) the multiple of the two , of an epicenter block in the nucleation process in the weak frictional instability regime. The model parameters are , , , , and . The time origin is set to the point where the epicenter begins to move. The arrow indicates the point of . The dotted horizontal lines represent the lines corresponding to , and to . The integers attached to the curves indicate the number of moving blocks . The insets are magnified views of the region around .

In order to see how the nucleation process evolves with the time, we show in Fig.6 typical time evolutions of (a) the block sliding velocity , (b) the state variable , and (c) the multiple of the two quantities , of an epicenter block in a typical nucleation process of a large event realized in the stationary state. The origin of the time () is set to be the onset of the nucleation process of the event. The model parameters are set to , , , and . The inequality is well satisfied, indicating that the system is in the weak frictional instability regime. The nucleation length estimated from eq.(30) to be given below is . The discreteness of the model is eminent in this regime.

At an early stage of the nucleation process, only an epicenter block moves. After some time, the neighboring blocks join this move one by one, causing a spatial expansion of the nucleus. As can be seen from Fig.6(a), the velocity of an epicenter block exhibits a step-like behavior, i.e., it exhibits an almost discontinuous rise when the block contingent to the moving blocks begins to move joining the nucleation process. As here, the system gets into the acceleration phase as soon as the number of blocks is increased from 3 to 4, and the epicenter-block sliding velocity begins to increase sharply. The block motion in the subsequent acceleration phase will be examined in the next subsection (Fig.6 to be continued to Fig.8).

One important general observation is that the epicenter-block sliding velocity in the initial phase stays very low up to , of order the pulling speed of the plate . This property can also be derived analytically as shown in §IVA below. In real faults, the plate motion is extremely slow, a few [cm/year] 1 [nm/sec]. Real-time detection of such a slow sliding motion would practically be impossible.

The state variable of an epicenter block initially takes a large value as shown in Fig.6(b). This is simply because linearly increases during the interseismic period according to eq.(7), acquiring a large value just before the onset of the nucleation process. During the initial phase, still keeps its large value since the velocity is still small on the r.h.s. of eq.(7), while it drops steeply beyond . The quantity increases with the time beyond , as can be seen from Fig.6(c).

We note that the dynamics of the model as shown here does not change much depending on the -value or on the -value as long as is taken smaller than , though the time evolution tends to be milder for smaller or larger . This tendency can naturally be understood because the smaller or the larger in eq.(6) means a larger contribution of the velocity-strengthening -term. The velocity-strengthening force serves to soften an abrupt change, causing a smoother time-evolution of observables.

The parameter choice of Fig.6 corresponds to and the discreteness of the model tends to be important around . In order to examine the effect of the discreteness on the nucleation dynamics, and to examine an approach to the continuum limit, we show in Fig.7 the corresponding figures for a different set of the parameters corresponding to Figs.4 and 5, i.e., , , (other parameters are the same as in Fig.6), which yields a larger -value of . One can see, while the nucleation process becomes smoother in this near-continuum case as expected, qualitative features remain similar to those observed in the strongly discrete case of Fig.6.

Figure 7: The time evolutions of (a) the sliding velocity , (b) the state variable , and (c) the multiple of the two (c), of an epicenter block in the nucleation process in the weak frictional instability regime. The model parameters are , , , , and , corresponding to the near-continuum regime. The time origin is taken at the point where the epicenter begins to move. The arrow indicates the point of . The dotted horizontal lines represent the lines corresponding to , and to . The insets are magnified views of the region around .

A characteristic feature of the block motion in the quasi-static initial phase is that there exist two different time scales: a slow motion of the time scale and a faster one of the time scale . The former might be better described by the slow time variable . Indeed, a perturbative treatment to be given in §IVA yields the time evolutions of the sliding velocity and of the state variable of the epicenter block as

(11)
(12)

where are constants to be determined by initial conditions, is the value of , and

(13)

with defined by

(14)

In the solution, the number of simultaneously moving blocks (the nucleus size) is assumed to be fixed during the block movement.

When the number of moving blocks or the nucleus size is small such that , both and are negative. When the condition is reached, changes its sign, leading to the instability. In fact, this condition determines the point of . From eq.(11), one can show that the block sliding velocity stays of order throughout the initial phase up to .

iii.3 C. The acceleration phase of the weak frictional instability regime

Next, we proceed to the acceleration phase of the weak frictional instability regime, which occurs beyond succeeding the initial phase. In the acceleration phase, the block movement exhibits a prominent acceleration, no longer quasi-static nor reversible.

In Figs. 8 and 9, we show typical time evolutions of (a) the sliding velocity , (b) the state variable , and (c) the multiple of the two , of an epicenter block. The model parameters are taken to be the same as those of Figs.6 and 7, the former corresponding to the strongly discrete case, and Fig.7 to the near-continuum case. The origin of the time () is set here to the point of . Note the difference in the time scales from those in Figs. 6 and 7: The abscissa in Figs.8 and 9 is , instead of in Figs.6 and 7. The arrows in the figures indicate the points of and of , respectively.

Figure 8: The time evolutions of (a) the sliding velocity , (b) the state variable , and (c) the multiple of the two , of an epicenter block in the acceleration phase in the weak frictional instability regime. The model parameters are , , , , and , corresponding to the strongly discrete case. The time origin is taken at the point of . The arrows indicate the points of and of . The dotted horizontal lines represent the lines corresponding to , and . The blue line in (a) is the theoretical curve for fixed , eq.(35).

As can be seen from Figs. 8(a) and 9(a), the epicenter-block sliding velocity increases rapidly in the acceleration phase, reaching a maximum of order , then decreases sharply and finally stops around . The state variable, which stayed nearly constant keeping its large value of order throughout the initial phase, begins to drop in the acceleration phase, and eventually becomes of order unity. Since the increase in dominates over the decrease in at an earlier stage of the acceleration phase, increases for some period, reaches a maximum, then drops sharply until it becomes close to unity: See Figs.8(c) and 9(c). Note that, around , the time variation of tends to level off exhibiting a much slower time dependence as can be seen from the inset. It is an inevitable consequence of the equation of motion, eq.(7).

Figure 9: The time evolutions of (a) the sliding velocity , (b) the state variable , and (c) the multiple of the two , of an epicenter block in the acceleration phase in the weak frictional instability regime. The model parameters are , , , , and , corresponding to the near-continuum regime. The time origin is taken at the point of . The arrows indicate the points of and of . The solid curves in (a) are the theoretical fitting curves, eq.(42) with eq.(44) (), and eq.(55) with eq.(57) ().

To a good precision, the maximum sliding velocity is reached when the relation is met. Eq.(7) indicates that, when the condition is met, the state variable takes a minimum. Meanwhile, sticks to a value close to unity in this range, yielding the relation . It means that the sliding velocity takes a maximum at the point where the condition is met.

In Fig.10(a), we show the time evolution of the number of simultaneously moving blocks, i.e., the nucleus size . The data exhibit a sharp peak at which the number of simultaneously moving blocks becomes maximum. This point was taken in subsection A as our tentative criterion of (). At or very close to this point, the epicenter block ceases to move (the double arrow in the figure), beyond which the group of simultaneously moving blocks are detached into two parts, each part propagating in opposite directions.

The point where takes a value unity and the epicenter-block sliding velocity reaches its maximum, might also be used as a reasonable criterion of . This definition of tends to yield a -value somewhat smaller than our previous definition of (), i.e., the maximum of the number of the simultaneously moving blocks. One justification of the new criterion might be the observation that the epicenter-block motion in the time range after levels off around has already become similar to the one observed in a typical block motion in the high-speed rupture phase. In this sense, the high-speed rupture has already set in in the epicenter region when the epicenter blocks satisfies the relation . Thus, in the following, we adopt as our criterion of the relation being reached at the epicenter block. This point agrees with the point of taking a minimum, or taking a maximum. In fact, the -values indicated in Figs.1,4,5,8 and 9 above were the ones defined in this way unless otherwise stated.

In earthquake dynamics, there generally exist two different types of velocities. One is the fault sliding velocity (particle velocity), corresponding in our model to the block sliding velocity . The other is the rupture-propagation velocity (phase velocity), corresponding in our model to the propagation speed of the rim of the rupture zone . In the nucleation process, the latter also coincides with the growth speed of the nucleus size ( half of it). Although the definition of the rupture-propagation velocity is somewhat obscure in the discrete BK model especially in the strongly discrete case, it might be well-defined in the near-continuum case as a (coarse-grained) growth rate of the rim of the nucleus. Namely, if the rim of the nucleus moves from the block to in a unit time interval, the rupture-propagation velocity might be defined by . We show in Fig.10(b) the time evolution of the rupture-propagation velocity computed in this way in the near-continuum case.

In the acceleration phase between and , we identify two characteristic points where the block motion appears to change its behavior. One is the point where the epicenter-block sliding velocity exceeds the crossover velocity , across which the -term gradually changes its character. The other is the crossover velocity at which the inertia effect becomes important. The inertia effect as meant here is borne by the first term of the r.h.s. of the equation of motion (6) or (8). This term tends to suppress the rapid acceleration, giving rise to the saturation and the subsequent drop of the sliding velocity . These two characteristic points also manifest themselves in our theoretical analysis of §IV below.

Figure 10: The time variations of (a) the total number of simultaneously moving blocks (the nucleus size) , and of (b) the rupture propagation velocity . The model parameters are , , , , and corresponding to the near-continuum case (the same as in Figs.7 and 9). The arrow indicates the point of , while the double arrow indicates the point where an epicenter block stops. The dotted horizontal lines in (b) represents the lines and . The solid curves in (b) are the theoretical fitting curves, eq.(46) with eq.(44) () and eq.(59) with eq.(57) ().

One sees from Fig.10(b) that the rupture-propagation velocity grows exponentially with the time until around , beyond which it grows faster than exponential (super-exponential). By contrast, as can be seen from Fig.9(a), the epicenter-block sliding velocity exhibits a faster-than-exponential growth even in the acceleration phase at . Namely, the sliding-velocity accerelation dominates over the nucleation-size expansion. Meanwhile, the super-exponential rapid growth of both the sliding velocity and the rupture-propagation velocity tends to be suppressed beyond the crossover velocity , which is caused by the inertia effect borne by the first term of eq.(8).

In Fig.11, we show (a) the epicenter-block sliding velocity , and (b) the rupture propagation velocity , versus the nuclear size normalized by , , instead of the time . The theoretical curves to be derived in §IV are also shown in the figure for comparison. From this, the changes of behavior at and at are eminent. The comparison with the analytical results are sometimes more direct in this form.

Figure 11: (a) The epicenter-block sliding velocity , and (b) the rupture-propagation velocity , plotted versus the nucleus size divided by the nucleation length, , in the acceleration phase in the weak frictional instability regime. The model parameters are , , , , and , corresponding to the near-continuum case. The time origin is taken at the point of . The dotted horizontal lines represent the lines corresponding to and . The theoretical fitting curves are also shown, (a) eq.(42) () and eq.(55) (), and (b) eq.(46) () and eq.(59) ().

iii.4 D. The acceleration phase of the strong frictional instability regime

Next, we study how the dynamics evolves during the acceleration phase for the case of the strong frictional instability. Remember that the model in the strong frictional regime lacks in the quasi-static initial phase.

The block motion here turns out to be similar to that of the weak frictional instability regime with a stronger discreteness. In Fig.12, we show the time evolutions of (a) the sliding velocity , (b) the state variable , and (c) the multiple of the two , of an epicenter block in the acceleration phase in the strong frictional instability regime. The parameters are taken to be , , , , and . In the event of Fig.12, the nucleation length is , i.e., the condition has been met during the one-block motion. Note that, in contrast to the weak-frictional instability case, this one-block motion is already irreversible.

One noticeable feature appears at an earlier stage of the high-speed rupture phase in the strong frictional instability regime. Namely, the block velocity often exhibits prominent oscillations with the time. The maximum sliding velocity realized at each oscillation is pretty high, comparable to that of a mainshock. In the inset of Fig.12, we show an expanded view around . Such an oscillatory behavior is rarely seen in the case of the weak frictional instability.

A closer look of the color plot in the inset of Fig.3 might reveal that such a velocity oscillation of the block velocity is borne by the propagation and the multiple reflections of the rupture front originally ejected at from the epicenter block. This rupture front propagates along the fault with an elastic-wave velocity , eventually becomes a rupture front of a mainshock. In the early stage of the high-speed rupture phase, this propagating rupture front is reflected every time it reaches a neighboring block, generating the second, third, rupture fronts, forming an oscillatory pattern. The period of oscillation should be given by , which, in the example of Fig.12, yields 0.5. This period is expected to be independent of the parameters like the plate loading velocity or the crossover velocity , while it might increase weakly with the frictional parameter because of the slowing-down effect due to the friction. Because of the friction, the velocity of the subsequent rupture-front propagation tends to be reduced, making the oscillation period a bit longer at a later time. We find that such expectations are consistent with the observation. For example, in an example shown in Fig.12, the observed oscillation period is 0.8-1.5, a bit longer than the expected value of 0.5.

Thus, in the strong frictional instability regime, the beginning of the high-speed rupture phase seems to be characterized by the multiple reflections of the propagating rupture front originally ejected from the epicenter site. After some time, the leading propagating rupture alone survives and propagates with an elastic-wave velocity for the major part of the mainshock.

Figure 12: The time evolutions of (a) the sliding velocity , (b) the state variable , and (c) the multiple of the two , of an epicenter block in the acceleration phase in the strong frictional instability regime. The model parameters are , , , , and . The time origin is taken at the point where the epicenter begins to move. The arrow indicates the point of . The dotted horizontal lines represent the lines corresponding to , and . The solid line in (a) is the theoretical curve for fixed , eq.(35). The insets are magnified views of the region around .

Iv IV. Analytical treatments

In this section, we wish to report on the results of our analytical treatments of the dynamical properties of the model, those based on either the perturbation theory [A]-[D], or on the mechanical stability analysis [E]. Readers interested only in the simulation results might skip to §V.

iv.1 A. Perturbation theory

We begin with the equation of motion (8) for the velocity variable . As mentioned, the low plate pulling speed provides an extremely small number, say, . Furthermore, throughout the nucleation process, the state variable tends to keep a very large value of order . Then, it might be convenient to introduce the reduced state variable of order unity, , defined by . One gets a set of equations of motions in terms of and ,

(15)
(16)

Then, we introduce the “first Fourier-mode approximation”, which states that for the most part of the nucleation process the spatial form of observables, i.e., the -dependence of the block sliding velocity or the block displacement , is given by that of the first Fourier mode. Namely, when the total blocks from to are moving, or is proportional to . In this approximation, it is implicitly assumed that the nucleus keeps a highly symmetrical form and the central block is an epicenter block.

We show in Fig.13 the spatial form of the block sliding velocity observed in our numerical simulations at several representative points of a typical nucleation process in the weak frictional instability regime, including (a) the initial phase, (b) the acceleration phase at , and (c) the acceleration phase at , together with the first Fourier-mode forms. Except for the time range beyond close to of Fig.(c), this approximation turns out to be reasonably good, allowing one to reproduce the motion of an arbitrary block within the nucleus by tracing only the motion of the central block (for odd ).

Figure 13: The time evolutions of the spatial pattern of the block-sliding velocity within the nucleus, (a) in the initial phase, (b) in the acceleration phase at , and (c) in the acceleration phase at . The model parameters are , , , , and . The solid curves represent the ones expected in the first Fourier-mode approximation where the height (the maximum velocity) is adjusted to the observed value.

Under this first Fourier-mode approximation, the equations of motion for the central block is given by

(17)
(18)

where has been given in eq.(14), and the subscript is dropped here and below.

Now, we expand and with respect to a small quantity , i.e., we perform a perturbation expansion in ,

(19)
(20)

At the zeroth order in , one gets a set of equations

(21)
(22)

The equation (21) for the zeroth-order velocity has two types of solutions, i.e., [A] , and [B] . The solution [A] describes the situation where the block is at rest when the plate drive is tuned off (). By contrast, the solution [B] describes the situation where the block is moving even when the plate drive is turned off. Hence, one expects that the solution [A] describes the initial phase, while the solution [B] describes the acceleration phase. In the following subsections, we analyze each case separately in some more detail.

iv.2 B. The initial phase

Here the block sliding velocity at the zeroth order is zero, . As was seen in §III, an eminent feature of the block motion in the initial phase is that there exist two time scales: One is a slow motion of order the loading velocity , and the other is a faster one of order unity. One way to deal with such two different time scales within the perturbative scheme might be to introduce two kinds of time variables, associated with a slow motion and associated with a fast motion, and regard various observables as a function of both and like , and . The original time derivative in the equation of motion is replaced by . Since , the zeroth-order equation (22) is reduced to . This means that depends only on , not on , i.e., .

The equations of motion at then reads as

(23)
(24)

Since the zeroth-order quantity is bounded in the limit, the corresponding first-order quantity needs to remain finite in the limit in order that the perturbation analysis remains meaningful, i.e., the relation is required in the limit. This relation is met if the equality

(25)

holds. From eq.(23), is obtained as

(26)

Substituting this into eq.(24) and taking the limit, one gets an equation to determine the hitherto undetermined -dependence of as

(27)

This can be solved to yield,

(28)

where . Substituting this into eq.(23), one can get a full solution of eq.(23) as

(29)

where are numerical constants to be determined via the initial condition, and is given by eq.(13).

While is always negative, is either negative or positive depending on whether or . When , both terms vanish quickly in eq.(29), whereas, when , the term grows quickly leading to the instability. In fact, the borderline case represents the nucleation length discriminating the stable nucleation process corresponding to the initial phase and the unstable nucleation process corresponding to the acceleration phase.

Now, from the condition , we reach an analytical expression of as

(30)

In the discrete BK model, the initial phase realized at is ever possible only when is greater than the lattice spacing or the block size, i.e., . Then, the condition of this nucleation length being greater than the block size yields the condition of the weak frictional instability,

(31)

where the quasi-static nucleation process is realizable in the discrete BK model. In other words, when , is less than the block spacing and the quasi-static nucleation process cannot be realized in the BK model due to its intrinsic discreteness. This is exactly the point discussed by Rice Rice (). Hence, either the weak or the strong frictional instability is determined by the relation between the two parameters and only, a strong instability for and and a weak instability for .

We emphasize that the continuum limit of the model corresponds to so that the continuum limit of the BK model with spatially homogeneous parameters always lies in the weak frictional instability regime, which accompanies the quasi-static nucleation process. Another derivation of and based on the mechanical stability analysis will be given in the following subsection E.

We note in passing that the analytic formula of given by eq.(30) is in excellent agreement with the -value determined numerically by artificially sopping the external loading as explained in §IIIA. Precisely speaking, the -value determined by artificially sopping the external loading could slightly deviate from the analytical result. Two reasons of such a deviation are identified. In one, a nonzero loading speed sometimes causes an “overshooting” giving a bias toward the instability. In the other, the spatial pattern of the block displacement and the block sliding velocity within the nucleus sometimes deviate from the one assumed in deriving the analytic form of §IV, i.e., of the first Fourier-mode form.

iv.3 C. Acceleration phase at

Next, we perform the perturbation analysis of the acceleration phase. Here, the block motion is no longer stable nor quasi-static, but is essentially unstable and irreversible. There is no slow process so that no need to consider . In contrast to the initial phase, the zeroth-order velocity describing this regime should be nonzero (the solution [B]).

We divide our analysis of the acceleration phase into the two time regimes from the technical reason, i.e., the regime of and of . In this subsection C, we deal with the regime . The regime will be dealt with in the next subsection D. For , eq.(21) reduces to the linear differential equation of the form,

(32)

whose solution is given by

(33)

where has been given by eq.(13). Since in the acceleration phase, is positive leading to the instability. The time evolution of the state variable is given by

(34)

In the analysis, the size of the nucleus is assumed to be fixed. Of course, an important part of the nucleation process, particularly in the unstable acceleration phase, is how the nucleus size expands with the time and how various observables evolve under the spatial expansion of the nucleus. In order to deal with such a nucleus expansion, we need additional information about the condition under which the block contingent to the moving blocks located at the rim of the nucleus begins to move. This condition actually depends on the stress state of the block assembly at the beginning of the nucleation process in question, which was basically set by the previous large event preceding the event in question.

We find from our numerical simulations that, in the steady state of an earthquake sequence, the excess stress , which is defined as the elastic-force difference at a given block between the initial value at the beginning of the nucleation process and the threshold value at which that block eventually begins to move involved into the nucleation process, is more or less constant over blocks involved in a given event, even though this quantity is scattered considerably over various events in an event sequence. This feature originates from the fact the stress distribution after a large event tends to be flat over blocks involved in this event.

Equivalently, the threshold displacement , which is defined as the displacement that the block located at the rim of the nucleus exhibits in order for the neighboring block initially at rest begins to move, also turns out to be more or less constant over blocks. In fact, there is a relation . In Fig.14, we show typical distributions of divided by its average over blocks involved in a given event , , for various parameter sets. The data for each parameter set is an average over events. As can be seen from the figure, tends to obey a common distribution characterized by a single-peak structure, suggesting that the approximation to regard (or ) to be constant over blocks involved in an event may not be so bad.

Figure 14: The distribution of the threshold displacement divided by its average over blocks involved in a given event , , for various parameter sets given in the legend. The other parameters are fixed to and .

For convenience of the description, we introduce the reduced time variable for which the time origin is taken at the point where the -block movement begins. In the symmetric block motion of the first Fourier-mode type we are considering here, the two blocks contingent to the nucleus begin to move entering into the nucleus motion of the size at the reduced time . We consider the series of nuclear sizes . Let the sliding velocity and the displacement of the central block at the transition from the -block motion to the -block motion be and . From eq.(33), one has

(35)
(36)

Within the first Fourier-mode approximation, the displacement of the block located at the rim of the nucleus means the displacement of of the central block. Hence, the the -constant condition for the transition can be given as the condition for the central block,

(37)

which, together with eq.(36), yields the equation to determine

(38)

Eqs.(35) and (38) yield the recursion relation for ,

(39)

which is solved as

(40)

To proceed further, we consider the situation where is large enough, . For , , and