A possible rupture process of slow earthquakes on a frictional fault

Summary. A possible mechanism for the occurrence of slow earthquakes is investigated by calculating numerical solutions for the dynamical rupture process on a quasi-three-dimensional fault with heterogeneous frictional strengths. Experimental friction laws for the dependence of sliding frictional stress on slip velocity, which are based on the cohesive properties of fault asperities, are taken into considerations. It is found that the applied stress does not drop very rapidly with time and the rupture velocity remarkably decreases as the dependence on slip-velocity becomes smaller. These deceleration effects for the rupture propagation are greatly enhanced with increasing heterogeneities in the distribution of frictional strength and as the initial shear stress has lower levels with respect to the average strength. For these cases, the growth of rupture is extremely slow in a nucleus region with the dimension as large as 10 times the initial rupture length, and gains a terminal velocity dependent on the above factors. The displacement-time function becomes noticeably extended in these cases, and indicates a stick-slip-like phenomena in the extended time interval for a strongly heterogeneous fault. It seems that these results could explain the characteristic features of slow earthquakes.

A rupture process of slow earthquakes 131 static and sliding frictional stresses are non-uniformly distributed and being subjected to a finite shear stress. These non-uniformities have been introduced to represent heterogeneous properties of fault material which has asperities or high strength barriers on the fault surface. This model provided satisfactory explanations to the initiation, spreading and stopping of rupture, stick-slip instability, radiation of high-frequency seismic waves, and non-uniform fault displacements, as well as to the occurrence of a main shock, multiple events, aftershocks on a single fault, and its statistics including magnitude frequency relationship (Mikumo & Miyatake 1979, henceforth Paper 11), while it did not account for aseismic fault slip. In this model, sliding frictional stress was assumed to be a function of position on the fault, but independent of time, slip displacement and slip velocity. If we introduce, however, non-elastic fault material such as a viscous layer or viscoelastic zone into the fault, the stress should be dependent on slip displacement and/or velocity. Dynamical solutions for these cases can be obtained in a similar way to that in Paper I, but it seems rather difficult to estimate the relevant non-elastic parameters from laboratory experiments.
The approach taken here is based on the characteristic friction laws on rock materials. It has been suggested (Weertman 1969) on the basis of friction experiments on metals (Bowden & Tabor 1964) that the frictional stress required to produce slippage across the fault interface decreases with increasing slip velocity. The velocity dependence of sliding friction is also found in friction experiments on silicate rocks (Scholz & Engelder 1976;Dieterich 1978). These friction characteristics are explained (Scholz & Engelder 1976) by the role of indentation and ploughmg of asperities on the sliding surface; the process that yields stick-slip instability in the frictional sliding is attributed to creep of indentation asperities on the surface, resulting from stress corrosion cracking or hydrolytic weakening, whch leads to the increase in real area of contact with time; t h s in turn yields increases in friction with increasing time of stationary contact or decreasing sliding velocity (Scholz & Engelder 1976). This feature is consistent with a slightly modified form of the Bowden & Tabor (1964) adhesion theory. From the laboratory data shown in Fig. 1, the following empirical forms have been given to relate the coefficient of sliding friction Pd to slip velocity V , (1 1 p d = & + f f log(bDc/T/+ 1) (Dieterich 1978) ( 2 ) (Scholz & Engelder 1976) analogous to the time dependence of static friction, although a more complex form has been proposed later (Dieterich 1979) to explain details of different experiments. po and P d are the dynamic friction coefficients at a specific sliding velocity and at the average velocity, respectively, and a, /3 and b are numerical constants. For several experiments, these constants have been estimated as; / 3 = 0.07 f 0.02, P d = 0.49 f 0.03, a! = 0.02 and b = 1 .O-2.0. D, is defined as the critical dimension of the contact of asperities and hence D,/V is regarded as an average life time of a population of contacts (Dieterich 1978). These two forms are essentially similar except for higher slip velocities. The experimental results of Scholz & Engelder (1976) shown in Fig. 1 indicate weak displacement dependence of sliding friction for smaller displacements in the case of constant velocity. Since the sliding friction should level down from static friction, we tentatively use a slightly modified form, introducing slip displacement D instead of Dc, where k estimated from the laboratory data is of the order of 0.2-0.4. The sliding frictional stress c7d in this case may be written, by multiplying a normal stress on to equation (3), as, where od = Eldon, us = psun and B = 00,. The variation of Ud with time is, -=B (x dt and its incremental increase can be calculated by,

T. Mikumo
where f(D, V ) simply replaces the right side of equation ( reason, we are dealing only with the velocity dependence in the subsequent numerical solutions. od also varies depending on the location of the fault if static frictional strength is non-uniformly distributed. Now we introduce the above condition into the previous fault model in Paper I. The fracture criterion used here is a finite-stress criterion as in the foregoing model; once the initial shear stress exceeds the static frictional strength at each point on the fault, slip motion begins. In this case, the slip motion is resisted by the sliding frictional stress

134
T. Mikumo dependent on slip velocity I/, which is given by equation (4) with k = 0. The external shear stress is assumed to increase very slowly (by about 4 bar yr-') due to the increase in tectonic loading resulting from relative plate movements with a velocity of 10 cm yr-' (Paper 11). This increase, however, has little effect on the rupture propagation within the time interval now in consideration. Under these conditions, numerical solutions are calculated for dynarnical rupture propagation in various possible cases. The parameters specifying the fault model and other constants are taken as almost the same as in Paper I, except that Az = 8 km to give larger displacements. The adopted static and lowest sliding frictional stresses and the initial stress are tabulated in Table 1.
3 Numerical solutions for slow rupture process In this section, several representative results of numerical calculations are presented, which include the time-distance relation and the pattern of rupture propagation, and the displacement time function, all of which indicate the natures of slow earthquakes.

C A S E O F H O M O G E N E O U S S T A T I C F R I C T I O N A L S T R E N G T H S
We have worked out 30 cases. In these cases, the static frictional strength is assumed to be constant over the fault plane except at the starting point of rupture. The assumed strengths are 202, 210, 230, 250 and 300 bar, and the lowest sliding frictional stress is taken as 120 bar. The initial shear stress i s also assumed to be homogeneous, being 200 bar or slightly larger. The B-value indicating the velocity dependence of the sliding frictional stress is taken here as ranging from 4 to 100 bar. When uo is close to us and also B is large, the case simply reduces to that treated in Paper I.
In the normal rupture process, the rupture front propagates with an elliptic shape with a nearly P-wave velocity in the direction parallel to the applied shear stress, and with a nearly S-wave velocity in the direction perpendicular to it, as shown in Paper I. For the fault dimension of 40 x 40 km, the rupture initiating from the centre of the fault reaches the prescribed fault edges in 3.7 and 6.3 s in the two directions, respectively. To compare with this normal case, an example of the rupture pattern in the present calculations is given in Fig. 3. Numerals on the rupture fronts indicate the propagation time elapsed after the rupture initiation. We see that in Case F with us = 250 bar it takes a much longer time for the rupture to cover the entire fault surface as compared with the normal process. Also noticeable here are the effects of decreasing B-values, which indicate the increase in the rupture propagation time from about 12 to 60 s at the left end of the fault as B decreases from 50 to 4 bar. The effects will be discussed later in more detail. Fig. 4(a)-(d) gives the time taken for rupture to propagate to different distances, for four different values of us with a parameter of B. The x-and y-axes here are taken parallel and perpendicular to the applied shear stress and measured from the starting point of rupture. From these figures, the following obvious features are immediately noticed. For fixed us and B-values the rupture propagation takes more time in the y-direction than in the x-direction, indicating an elliptical mode of propagation as expected from the previous results. It has been pointed out, however, that this is not always the case in the normal rupture process on a more complete three-dimensional fault (Miyatake 1980); if the static frictional strength is considerably larger than the applied stress, the rupture tends to show a circular propagation (Miyatake 1980) with decreased velocities in the x-direction. For this reason, only the curves in the y-direction is presented for a, = 300 bar. * .
I .

.
* . It is also noticed from a comparison between Fig. 4(a)-(d) that for a fixed B-value the rupture propagates more slowly as the static frictional strength us has higher levels with respect to the applied stress uo. When us is elevated from 210 to 300 bar, the propagation time of rupture to reach a position at y = 20 km increases from 5.6 to 18.6 s for B = 50 bar and from 15.6 to 97.5 s for B = 4 bar. This type of variation in the rupture velocity has been suggested analytically in two-dimensional shear cracks (e.g. Burridge 1973; Andrews 1976) and numerically in a three-dimensional fault (Miyatake 1980), although these studies do not involve the effect of velocity dependence of sliding frictional stress. Although the above decrease in the rupture velocity is not a new finding here, the difference in the stress level could be one of possible sources for slow earthquakes, and the effect would be enhanced if combined with other sources.
The most remarkable feature to be noted here is that for a fixed value of us the propagation time of rupture rapidly increases with decreasing B-values. The rates increase by a factor of 2.8, 4.8 and 5.3 respectively in the case of us = 210, 250 and 300 bar when B decreases from 50 to 4 bar. This clearly indicates that small B-values greatly decelerate the rupture propagation. This effect results from the assumed cohesive properties of fault material due to ploughmg asperities on the fault surface. From the laboratory data for a  limited number of rocks (Scholz & Engelder 1976), the B-value in this case may be of the order of 5-20 bar for us = 200-300 bar, which is small enough to cause slow rupture propagations. Another important feature is the transient nature of rupture velocities. As is clear from Fig. 4(a)-(d), the rupture moves very slowly at the initial stage, then gradually increases its propagation velocities up to a distance of about 6-10km, and finally gains a terminal velocity, for all cases of us and B-values. It has been stated (Andrews 1976) that a planestrain shear crack with a slip-weakening law shows slow initial growth of rupture up to a distance comparable to the doubled initial crack length, with its terminal velocity approaching the Rayleigh velocity. To compare with this, the initial and terminal velocities in the y-direction perpendicular to the applied stress are explicitly shown in Fig. 5, as a
-. ^^ . -..  ^I. --. -" --. ---. " ^-I"^ .  in this case does not seem to be dependent on the stress level, but might be regarded as the nucleus of a slow earthquake, where strain energies are accumulated for further accelerating ruptures. It is not clear, at this moment, what conditions specify the nucleus area. delayed arrival of the rupture to the centre of the fault (x = 20 km,y = 20 km). It should be remarked that the rise time at this point appears considerably shorter than that at the rupture initiating point unlike in the normal rupture process.

E F F E C T O F H E T E R E O G E N E O U S D I S T R I B U T I O N O F S T A T I C A N D S L I D I N G F R I C T I O N A L S T R E S S E S
We have made a number of calculations for the various cases given in Table 1, in order to examine the effects of spatially heterogeneous distribution of static and sliding frictional stresses. Three different distributions are shown in Fig. 7, which are similar to those tested in Papers I and I1 (Mikumo & Miyatake 1978. Case A includes a number of high strengths exceeding 500 bar within the encircled portions, Case B has smaller inclusions with higher strengths, and Case C has a weakly heterogeneous distribution with a standard deviation of 25 bar. The average strength in these cases is 3 19,282 and 260 bar, respectively. The sliding frictional stress in these cases is also non-uniformly distributed as indicated by equation (4), and their minimum value is taken as 100 bar in Case A and 120 bar in Cases B and C. The initial stress and B-values are taken as the same in the foregoing section.  Fig. 3. A close comparison between the two cases for B = 50 and 20 bar shows, however, that it requires a somewhat longer time in Case C than in Case F for the rupture to reach the left edge of the fault. For more heterogeneous cases, Cases A and B, the shape of the rupture fronts becomes strikingly complicated and its propagation velocity remarkably decreases, as high strength inclusion increases and as the B-value decreases. The rupture velocities in these two cases are given in Fig. 9, analogous to Fig. 5, where their possible ranges are shown by shaded zones and Vx(I), V,,(I) and V y ( T ) denote the average initial velocities in the x-and y-directions over a distance of l O k m and the terminal velocity in the y-direction, respectively. The decreases in the rupture velocity in these cases are again remarkable. In order to see more explicitly the effect of heterogeneities on the rupture velocity, Cases A and F with us = 300 bar, which have nearly the same level of the average strength, are compared in Fig. 10 for 8 = 6 bar. The most obvious effects found in Case A, as compared with the pattern in Case F, are the irregularly deformed rupture fronts and much longer propagation time of about 240 s to reach the left-bottom of the fault. It is thus demonstrated that strongly heterogeneous distribution of static and sliding frictional stresses also gives a dominant deceleration effect for rupture propagation. It can also be recognized in Case A that there is a nucleus zone where the initial rupture grows up in extremely slow speeds, although the initial velocity and the finiteness of the zone cannot be estimated correctly because of their appreciable difference in different directions. The displacement-time function in Cases A, B and C are shown in Fig. 11 (a), (b) and (c), where the location of the selected points are the same as in Case F. It is found that elongations in the total time duration and gradual decreases in its initial slope with decreasing B-values are extremely noticeable. The rise times of the function in the three cases are also depicted in Fig. 12, which clearly shows a rapid increase for smaller B-values indicating the steepest rate for Case A. If we compare Case A with Case F for us = 300 bar, and Case C with Case F for us = 250 bar, it is again noticed that the effects of heterogeneities appears remarkable in the displacement-time function. It has been shown in Paper I that the normal dynamic rupture process on a heavily nonuniform frictional fault yields stick-slip-like phenomena, which are manifested by episodically rising forms of the displacement time function in several seconds. For a most hetereogeneous fault, Case A in the present case, these minor fluctuations within a short time range are simply masked in the case of larger B-values, whereas for smaller B-values the stick-slip appears more elongated in a longer time range. Fig. 13 gives the displacement-time function in Case A for B = 6 bar. The rise time reaches 200 s at the initiating point of rupture and its vicinity. It is interesting to note that the displacements and stress drops in the nucleus region remain at low levels up to about 130-140 s and then start to develop into larger values. The earlier part might be regarded as a precursory slow movement to the later stage, although the entire displacement time function also shows characteristic features of a slow earthquake. The time which appears to separate the entire function into two parts corresponds approximately to the time by which the slow rupture reaches an edge of the prescribed fault boundaries. if this is true, it seems that propagating energies of the slow moving rupture are blocked by the boundaries and then turned back into the inside of the fault.

Discussion
It has been shown in the present paper that the velocity dependence of sliding frictional stress based on the cohesive properties of fault asperities, heterogeneous fault strengths, as well as the initial stress level with respect to the average frictional strength, would play a dominant role in generating slow earthquakes. On the other hand, Yamashita (1980) concluded that a low initial stress combined with low decreating rates in the difference between the initial stress and sliding frictional stress with respect to the distance from the rupture initiating point, and high specific fracture energy will yield slow earthquakes. In the present calculations, we have dealt with the case corresponding to the lowest rate with a constant initial stress, including other possible sources. The fracture energy for creating a new crack surface at the advancing crack tip may be related to the cohesive strength of fault material, and would have some relation with the present case.
Non-elastic properties of fault zone material suggested in several researches might also be able to provide some explanations. For the case of viscous friction on the fault surface, the sliding frictional stress ud would be Ud = uo -qV/d, taking into account the constitutive equation u = 76, where r ) is the viscosity, d is the thickness of a gouge layer, and i is the strain rate. When uo reaches the level of us, ud drops from us in a linear proportion to slip velocity V . If this is tentatively compared with the experimental relation (4) in the present case, q/d would correspond to B log V/V, and hence q would be of the order of 108poise for B = 10 bar, Y = 2-10 cm s-l and d = 10 m. This does not seem to be an acceptable value, although this oversimplified assumption does not necessarily rule out the possibility of non-elastic models.
It has been reported in the literature (Ando 1975) that the duration time of faulting for earthquakes recurrent on the same fault plane is not always constant as manifested in the 1854 Ansei I1 and 1946 Nankaido earthquakes. If this is true, a single fault could generate a slow earthquake in one period and a normal brittle earthquake in another. This feature might be explained by a temporal difference in the state of fault zone material such as partial melting and adherence. The present results may provide another interpretation; faulting motion during a large-scale earthquake would break the contact of asperities on the fault surface, and consequently the cohesive strength of plouglung asperities and heterogeneous frictional strength would become considerably different at the time of a succeeding earthquake when the tectonic shear stress again reaches the level of the lowest frictional strength.
The present numerical results show that unusually slow ruptures with small displacements and stress drops could grow up in a limited nucleus region on a strongly heterogeneous fault with smaller B-values and high average frictional strength with respect to the initial stress. The slow rupture develops into subsequent larger slip movements, indicating a stickslip-like phenomena, after the time it reaches the fault edges. These features appear consistent with those for laboratory friction experiments and their simulation including preseismic stable slip and transient unstable slip motion (Dieterich 1979). The time duration in the present study had to be limited to within 4 min, but if dynamic or quasi-static solutions could be extended to much longer times, the natures of preseismic slip and its triggering mechanism for a normal earthquake would be more explicitly elucidated.

Conclusions
We have investigated a possible mechanism to generate slow earthquakes through the dynamical rupture process on a quasi-three-dimensional fault with heterogeneous frictional strength, taking into consideration experimental friction laws for the cohesive properties of indentation and ploughing fault asperities. The effective factors to yield slow earthquakes we have found here are the dependence of sliding frictional stress on slip velocity based on the cohesive properties, which is indicated by a coefficient B, spatially heterogeneous distribution of static and sliding frictional stresses, us and ud, and the difference in the level between the applied shear stress uo and the average frictional strength us. The main results from the present numerical calculations are as follows: (1) For smaller B-values the sliding frictional stress and hence the initial stress does not drop very rapidly with time. For a fixed value of us, the rupture velocity remarkably decreases with decreasing B-values, which give appreciable deceleration effects.
(2) For a fixed B-value, the rupture velocity decreases as the difference in the stress level us-uo becomes larger.
(3) For smaller B-values and larger difference in the stress level us-uo, the growth of rupture is extremely slow in a region with a dimension of about six to ten times the initial rupture length. The rupture velocity then slowly increases and reaches a terminal velocity dependent of the B and us-uo values.
(4) Heterogeneous distribution of static and sliding frictional stresses has also noticeable deceleration effects on rupture propagation. As high strength inclusion increases, the rupture velocity goes down to extremely low values in a nucleus region.

T. Mikumo
(5) The displacement-time function becomes greatly elongated with decreasing Band increasing us-cro values and increasing heterogeneities. For a strongly heterogeneous fault, the function shows a stick-slip-like phenomena in an extended time interval.