22 August 1994 PHYSICS LETTERS A
ELSEVIER
Physics Letters A 191 (1994) 398-402
Scaling of synthetic seismicity from a one-dimensional dissipative, dynamic lattice model Jeen-Hwa Wang Institute of Earth Sciences, Academia Sinica, P.O. Box 1-55, Nankang, Taipei, Taiwan, ROC
Received 27 July 1994; revised manuscript received25 March 1994; accepted for publication 9 June 1994 Communicatedby A.R. Bishop
Abstract
Seismicity is synthesized from the one-dimensional dynamic lattice model of Burridge and Knopoff [ Bull. Seism. Soc. Am. 57 (1967) 341 ] with a fractal distribution of breaking strengths. A linearly velocity-dependent friction controls the sliding of the lattice point. The stiffness ratio (s) is a parameter to represent the level of conservation of the system. In a range ofs values, the relation of frequency versus magnitude follows a power-law function, and the scaling exponent (b) of the relation relates to s in the form b~s -2/3.
The relation of frequency (N) versus magnitude (M) of earthquakes displays a power-law function [ 1 ]. Numerous studies based on the quasi-static equivalent of the BK model (i.e. the lattice model proposed by Ref. [ 2 ] ) through the cellular automaton method were carried out for studying this scaling. Some authors (see Ref. [ 3 ] ) only considered conservative systems and ignored the dissipation effect. However, Nakanishi studied the scaling based on a 1-D model including a dissipation term [ 4 ]. His resuits showed that the scaling exponent increases somewhat with the conservation level of the system. The dynamic model must be more appropriate for the study of mechanical processes including dissipation of energy than the quasi-static one. For the 1-D homogeneous Burridge-Knopoff model with velocity-weakening friction considered in Refs. [ 5,6 ], and the cellular automata analog considered in Ref. [ 4 ], the scaling exponent in the magnitude versus frequency distribution was nearly constant, approximately - 1, as a stiffness parameter, representing the ratio of compressional to shear elasticity, was varied.
In this paper, an attempt is made to examine how the scaling law varies in a spatially inhomogeneous model as the analogous parameter is varied. The I-D BK model with a fractal distribution of breaking strengths is taken into account. The dynamic friction which increases with velocity controls the motion of the mass. Most strikingly, in this case it is found that the scaling exponent varies continuously with the stiffness ratio, defined below, more analogous to results obtained for a self-organizing cellular automata model with velocity-independent dynamic friction introduced by Ref. [ 7 ]. The BK model [2], which was originally applied to approach the fault dynamics, consists of a chain of N masses with equal mass ( m ) and springs with each mass being linked by a spring of strength (K) with two other neighbors and each mass also being pulled through a spring of strength (L) by a moving plate with a constant velocity (l/p). This system is illustrated schematically in Fig. 1. Initially, all masses rest in an equilibrium state and the spacing between two masses is a. Each mass is located at position un, mea-
0375-9601/94/$07.00 © 1994 ElsevierScience B.V. All rights reserved SSD10375-9601 (94) 00502-G
399
J-H. Wang/Physics Letters A 191 (1994) 398-402
".V
/ / / Fig. 1. O n ~ e n s i o ~ l
sured from its initial equilibrium position, along the x-axis. Each mass is subjected to a frictional force, F,. The equation of motion at the nth mass of the system is
miin=K(un+s-2Un +Un_l)-L(un-Vpt)-Fn.
(1) The dots represent differentiation with respect to time t. Although the spacing a is remarkably not an explicit parameter in Eq. ( 1 ), it is a length to provide the short wavelength cutoff. A parameter called the stiffness ratio (s) is defined as the ratio o f K t o L. In other words, s represents the ratio of the coupling between two masses and that between a mass and the moving plate. A larger s value denotes a stronger coupling between two masses than that between a mass and the moving plate and produces a smaller amount of loss of energy through the spring L, thus indicating a higher level of conservation. Therefore, s is a significant parameter to represent the level of conservation of the system. For simplification of computation, the m and L values are respectively taken to be unity, and thus the stiffness ratio s is simply equal to K. For the earth, Vp is very small, of the order of 10 -t2 units. But in practical computation, a value of 10 -4 units is taken to avoid the need for a long computational time. The boundary condition was considered by Ref. [ 7 ] to be a factor affecting the scaling exponent in the simulation. Nevertheless, only the periodic boundary condition is applied at the two end masses. The velocity-dependent friction was first found by Ref. [ 8 ]. Such a velocity-dependent friction includes two processes: ( 1 ) the velocity-weakening process as the sliding velocity is smaller than a critical velocity and (2) the velocity-hardening process as the sliding velocity is greater than the critical velocity. The critical velocity (re) is the one at which the dynamic fric-
m~u~-spring model.
tional force reaches the minimum value. T h e generalized velocity-dependent friction law is complicated. In this study, a linear friction law is considered and shown in Fig. 2, in which Fo is the static frictional force or breaking strength and gFo is the minimum dynamic frictional force. Fig. 2 shows that no backward motions are allowed or the sliding velocity is always greater than or equal to zero. This friction law shows a sudden drop of the frictional force from a static one to a dynamic one and the vc value is close to zero. Such a sudden drop of the frictional force behaves like an impulse to provide energy to the mass for further sliding. A smaller g value will produce a bigger force drop, thus it is easier to generate a largersize event. As the sliding velocity is greater than zero, the dynamic frictional force increases with sliding velocity, thus increasing the resistance to the motion. Hence, in addition to the coupling between the mass and the moving plate, the velocity-dependent friction must be also an important factor in the dissipation of energy of the system. A s t u d y of the friction effect will be presented in a separate paper. The fault zones where earthquakes occur are usu-
+F
+gF
+V©
Velocity
Fig. 2. Linearly velocity-dependent friction law:. Fo is the breaking strength; Vcthe critical velocity; g the friction drop ratio.
400
J-H. Wang/Physics Letters A 191 (1994) 398-402
ally quite complex. Map and field observations [ 9 ] and laboratory observations [ l0 ] showed that the distributions of physical and geometrical properties over the fault zone are usually fractal. Fractality is a very common property for natural phenomena [ 11 ]. Hence, the distribution of breaking strenghts of the system is considered to be inhomogeneous and follows a fractal function. The midpoint displacement method developed by Ref. [ 12 ] is used to yield a fractal distribution. From this method, only the fractal having Npoints ( N = 2 n+ ~) is generated at the nth computational level. In this study, since only the effect due to the change of the s value is considered, the g, n, and D values are taken to be constant and are 0.8, 7, and 1.5, respectively. The maximum value of the breaking strength used is 5 units. At a certain mass, as the sum of driving force due to the moving plate and spring forces from its neighbors exceeds the breaking strength, the mass is accelerated and starts to slide. Meanwhile, the breaking strength drops to the dynamic frictional force. After a while, the increase of both spring forces due to the change of the relative positions of the mass and its neighbors and the dynamic frictional force decelerates the sliding of the mass. Finally, the mass stops and sticks, and thus resulting in a drop of force. The moving plate will reload the mass to pull it to slide again. As a mass sticks after motion, the dynamic frictional force is considered to immediately recover to the static one, in other words, there is no delay time for the healing of friction. The displacement of a mass is measured from its new equilibrium position to the position where it sticks after motion. The position where the mass sticks after sliding is a new equilibrium position for the next stage of motion. Several connected masses may slide almost simultaneously during a time span. The sum of displacements of such masses in terms of time shows the time history of the displacement of such masses. Such a time history is considered to represent an event. The seismic energy (Es) is proportional to the maximum slip Um~, and the related force drop (Af), i.e. Es=um~Af. Therefore, the logarithmic value of the sum of the seismic energy of the masses for one event is taken to be the magnitude (M) of the event, i.e. M = l o g ( ~ Esi). Obviously, the present magnitude is an energy-based magnitude rather than the commonly used magnitude based on the peak amplitude of the seismogram.
However, for real earthquakes the former relates linearly to the latter. We integrated numerically Eq. (1) up to several 107 time units for numerous cases with different s values. The N versus M relation is derived from the synthetic events. As s < 10, the synthetic seismicity cannot show a good power-law distribution of frequency versus magnitude. As s> 120, only the delocalized events, for which all masses slide in a short time span, are generated. Hence, only the cases with s values from 20 to 120 are taken into account. However, the results for s = 5 and l0 are still given for comparison. Fig. 3 displays the spatial-temporal patterns of events for three models with s = 10, 30 and 50. The line segment linking up the slid masses represents an event. A longer line segment includes a larger number of masses and thus represents a larger event. For a smaller s value ( = 10 or 30), small and medium events from the majority of the synthetic events; while for a larger s value ( = 50), a greater number of larger events can be produced. It is evident that the number of larger events increases with s. Besides, the occurrence time of two events seems to be shorter at the masses with lower breaking strengths than at those with higher ones. The values of log N and M for 5 ~ Me, the log N value slightly increases with s. The log N values for 30 ~ s ~<100 vary only in a small range. This indicates that the change ofs value has different effects on the two regimes of small and large ones. A smaller s value makes the rupture be smaller, thus leading to more smaller events; while a larger s value makes the rupture be more capable of extending to a larger range, thus resulting in more larger events. The plot of s = 10 shows that only in a very small magni-
401
J-H. Wang/Physics Letters A 191 (1994) 398-402
k" ~" -,~.'f'~.' °," "-~".,-,. ,' :1"~,~. ';: o • ".~. ¢";.,~-"..: ~ :. "3 .'..~.. :1,~ .~. ~,- ~)i i'" f "' '~'""; " ;~.~::~"~'.'"~" .",'-" '~:"'"d"':~ ~.',~*"."~ :",:"~:, ",'-~':-' ~ "
.-,~:r. *i,'~ i
I',~ • ,e " . ' , ~ o " . '~ . v . ' , : . :~.. ~, ",~" ,~,~ 2.#,.-,.. -,-..'~ • ,, .-i', ".,':: . " ,'I: " ,:"." ~,t , ~ : ' : , ~ ' . " : - ' " - • , " : , ;" e ; t : ' $ : : • .~',: • ~ " " - , ' " I'.• 5 ,, ,, ~ . ' -.,,, .%...~,..,.~...-, .......... , .- './-- •. • -~,o. - . . '.",, , %'°, ~. .~¶~ , t ,',', ~,- , , - a I ~.'z ".." ~".) ..~.:'"...':".~-' "',-° " ~: ".."." "•';i ".)'~ "~. z: •d'-:s*o .'~;'";~" ",.'G'.tl ~t ~p ~l~ "f~.t• t : . ~ C o ; | i r o . ~ a- . f~,~Oo . a ' . $|.7. #)*.:" l e t o~ t. " . 6 " . ° o r.." oH b
[%[q.:;:.~,-., ~]...~.~-~,~, ,.:,.~,~ ,~.:~2.;-::.?.. <....¢.:.~'.~...,:~ ,-; ;.:,.,.:~ t '~/. ,e:. ~.,::t..~.$...~; ~1 .~;, , .5 e /,'¥, S ° • ~ %~- ~ f ~q" f%'.~ I* *" , " • ," i l ) , tI * ~.*, to { "
l~.:~*',.~,::-~'~.~',~;',.?,,...,,.,'. .;: ~" /-.;'o:.,~;~,Z.~'r ~'.•, ' v ' ; : . " ) - ~ )':.'1 ~..';~.Z,.';'::¢,,-,':
C
,,.-,~
,,,(.-x.
,~*~.~,:.?.':
.,.
a .&..,-,-
..,...t,:-:l
r "e:e{'~ .:~..,. ;~.::'~:..~ :j~ ]..!;~'~ :.';.~ ~:.-.# -.%~,t' ~: ",-~',.~ .~.,:r, ~ •;:,~Xo£ ~3 ":~,;~:," :J -
-t-~
.' ...
• ,~" • ~ :;.'..~ ". , . ,
'~-
~ •
.
~
~.
,~
[,t~.., ..,: -:,'.~,. J ,, -,,'.%:, .~ .' ,., ,'-'" , . ' . , ",.'.~'.. "~" ".,.:. t " .A., ~.. ~ t.':t ,';[| f.~",'. "r... i,1,. :l'~..t t ,t.~ ,,~::.;.' ,. ~.,. , - "., "..,. ',,.'.- . ' , .'," . ' . . ' ; . . " - - - .:-'1 l
J
II
l
o
t
|
|
7
units)
II
i
TIME
l
2
I (in
10
Fig. 3. The spatial-tcmpond patterns of ¢~mts: (a) for s = lO, (b) f o r s = 30 and (c) f o r s = 50.
4
0.5 0 0
(D Z
0 > _Q
d
0
-0.5 0
-3
Mognitude
3
0.5
2.5
s volue Fig. 5. b versus s in a logarithmic scale. The slope value of the solid line is about - 2/3.
Fig. 4. The data points of log N versus M: ( X ) s= 110, ( • ) s= 90, ( + ) s=70, (O) s=50, ([2) s=30 and (A) s=10. tude range the d a t a points follow a power-law function. O n the other hand, for s >i 20, the d a t a points o f log N versus M are d i s t r i b u t e d a l m o s t a r o u n d a line in a large m a g n i t u d e range a n d follow a power-law function. Fig. 4 also displays that the scaling exponent ( b ) , i.e. the slope o f the regression line o f the d a t a points, decreases with increasing s. In o r d e r to explore the relation between b a n d s, the plot o f b versus s for 5 ~_.20 are d i s t r i b u t e d almost along a line, whereas the d a t a points for s < 20 d e p a r t f r o m the linear t r e n d
very much. The regression line, inferred from the d a t a points for 20~
402
J-H. Wang~Physics Letters A 191 (1994) 398-402
nitude distribution functions were calculated by Ref. [4] for four cases with ,4=0.8, 0.85, 0.9 and 0.95 as s=2.00, 2.83, 4.50, and 9.50, respectively. These s values are smaller than the lower bound of s values used in this study. Although the b value was not given by Ref. [4 ], the results of Ref. [ 4 ] still showed that the scaling exponent increases with A or s. In Refs. [5,6], for three 12 ( = s ) values of 36, 64, and 144, which are almost within the range of s values used in this study, the scaling exponents estimated are approximately - 1 and not dependent on l 2 ( = s ) . In the cellular-automaton modeling by Ref. [ 7 ] for the 2-D models, for the nonconservative phase with a, = a2, the scaling exponent decreases continuously as a function of level of conservation a ( = or1= a2) for both free and open boundary conditions, ot equals s~ (4s+ 1 ) and increases with s. For both boundary conditions, as a < 0.24 the b value decreases almost linearly with increasing a; while as ot>0.24 and b value changes remarkably. The s values from 20 to 120 used in this study are nearly equivalent to the a values from 0.243 to 0.249 in Ref. [7 ]. Such a range of a values is just a very small part of ot values in Ref. [ 7 ]. Although only two b values were computed in such a range, the decreasing trend ofb with respect to s in Ref. [ 7 ] is similar to that in this study. As the a value is less than 0.05 for the free boundary condition and less than 0.01 for the open boundary condition in Ref. [ 7 ], the activity of the system does not follow a power-law function. The s values associated with a = 0 . 0 5 and a = 0 . 0 1 are about 0.0625 and 0.1667, respectively, which are much smaller than the low bound of the s value in this study. Evidently, the activities for the 2-D quasi-static BK model with a small s value, for which the coupling between two masses is much weaker than that between a mass and the moving plate, can show a power-law function. There are two main differences between this study and Ref. [ 7 ]. The first one is the way for simulation: a dynamic approach for the former and a quasi-static one for the latter, and the second one is that the dissipation term due to friction is included in this study but not in Ref. [ 7 ]. Although the dimensions of the models of this study and Ref. [ 7 ] are different, the comparison of the two results seems able to suggest two concluding points: ( 1 ) the loss of energy through
the L spring is a factor but not a main one for the nonconservation of the system to affect the scaling of the N versus M relation and (2) the dissipation of energy due to friction significantly reduces the level of conservation and rises the lower limit of the s value for this scaling. Hence, to retain the power-law scaling for the N versus M relation, a way to compensate the dissipation of energy due to friction is to strengthen the coupling between two masses. The author thanks the referee's useful comments. Thus study was financially supported by the Academia Sinica and the National Science Council (under Grant NSC 82-0202-M-001-154), Taipei, Taiwan, ROC.
References [ 1 ] B. Gutenberg and C.F. Richter, Ann. Geofis. 9 (1956) 1. [2]R. Burridge and L. Knopoff, Bull. Seism. Soc. Am. 57 (1967) 341, [3] P. Bak and K. Chen, Physica D 38 (1989) 5; P. Bak and C. Tang, J. Geophys. Res. 94 (1989) 15635; K. Chen, P. Bak and S.P. Obukhov, Phys. Rev. A 43 ( 1991 ) 625. [4] H. Nakanishi, Phys. Rev. A 41 (1990) 7086. [5] J.M. Carlson, J.S. Langer, B.E. Shaw and C. Tang, Phys. Rev. A 44 ( 1991 ) 884, [ 6 ] J.M. Carlson and J.S. Langer, Phys. Rev. A 40 (1989) 6470. [7] IC Christensen and Z. Olami, Phys. Rev. A 46 (1992) 1829; K. Christensen and Z. Olami, J. Geophys. Res. 97 (1992) 8729; Z. Olami, H.J.S. Feder and K. Christensen, Phys. Rev. Lett. 68 (1992) 1224. [8] J.H. Dieterich, J. Geophys. Res. 82 (1977) 5648; T. Shimamoto, Science 231 (1986) 711. [9] C.A. Aviles, C.H. Scholz and J. Boatwright, J. Geophys. Res. 92 (1987) 331; P.G. Okubo and K. Aki, J. Geophys. Res. 92 ( 1987 ) 345; C.H. Scholz and C.A. Aviles, in: Geophys. Mono. Ser. Vol. 37. Earthquake source mechanics, eds. S. Das, J. Boatright and C.H. Scholz (AGU, 1986) p. 147. [10] S.R. Brown and C.H. Scholz, J. Geophys. Res. 90 (1985) 12575; W.L. Power, T.E. Tunis, S.R. Brown, G.N. Boitnott and C.H. Scholz, Geophys. Res. Lett. 14 (1987) 29. [ 11 ] B.B. Mandelbrot, The fraetal geometry of nature (Freeman, San Francisco, 1983); D.L. Turcotte, PAEGEOPH 131 (1989) 171. [ 12] D. Sanpe, in: The science of fraetal images, eds. H.O. Peitgen and D. Sanpe (Springer, Berlin, 1988).