one-dimensional dynamical Modeling of Earthquakes : A review

Studies of the power-law relations of seismicity and earthquake source parameters based on the one-dimensional (1-D) Burridge-Knopoff’s (BK) dynamical lattice model, especially those studies conducted by Taiwan’s scientists, are reviewed in this article. In general, velocityand/or state-dependent friction is considered to control faulting. A uniform distribution of breaking strengths (i.e., the static friction strength) is taken into account in some studies, and inhomogeneous distributions in others. The scaling relations in these studies include: Omori’s law, the magnitude-frequency or energy-frequency relation, the relation between source duration time and seismic moment, the relation between rupture length and seismic moment, the frequency-length relation, and the source power spectra. The main parameters of the one-dimensional (1-D) Burridge-Knopoff’s (BK) dynamical lattice model include: the decreasing rate (r) of dynamic friction strength with sliding velocity; the type and degree of heterogeneous distribution of the breaking strengths, the stiffness ratio (i.e., the ratio between the stiffness of the coil spring connecting two mass elements and that of the leaf spring linking a mass element and the moving plate); the frictional drop ratio of the minimum dynamic friction strength to the breaking strength; and the maximum breaking strength. For some authors, the distribution of the breaking strengths was considered to be a fractal function. Hence, the fractal dimension of such a distribution is also a significant parameter. Comparison between observed scaling laws and simulation results shows that the 1-D BK dynamical lattice model acceptably approaches fault dynamics.


IntroductIon
Since the end of the 19 th century, several self-similar properties of spatial distributions and temporal variations of earthquakes have been found.These properties are generally characterized by a power-law function.Omori (1895) observed the temporal variation in the number n(t) of aftershocks following a mainshock in the form: n (t) t c p = , where c is a constant and p is the scaling exponent and is close to 1 for most observations.Utsu (1961) generalized Omori's law to a form of n (t) (1 t) c p = + .Gutenberg and Richter (1944) reported a frequency-magnitude (FM) scaling law of earthquakes in the form: logN = a -bM, where M is the earthquake magnitude and N is the cumulative or discrete frequency of events with magnitudes ≥ M. Ekstrom and Dziewonski, (1988) stated that the seismic energy, E s , released during an earthquake relates to M in the form: log(E s ) ~ ξM, where ξ is 1 for earthquakes with M s < 5.3 and 3/2 for those with M s > 6.8.Hence, there is a power-law function between N and E: N ~ E s -B , where B b p = .Kanamori and Anderson (1975) reported a power-law function for seismic moment, M o , and rupture length, L, of earthquakes in the form: M o ~ L 3 .This relation is commonly considered to hold up for a wide range of events (Hanks 1977).However, for Japanese earthquakes, Shimazaki (1986) reported a change in scaling from M o ~ L 3 to M o ~ L 2 occurring at the point where M o = 7.5 × 10 25 dyne-cm.Aki (1992) and Romanowicz (1992) also reported such a change for worldwide and California earthquakes.However, from more than 200 global earthquakes, Wang and Ou (1998) did not find such a change in scaling, and they also observed M o ~ L 2 in a large range of M o .
Source rupture duration, τ c , is an important parameters in the description of earthquake rupture (Kanamori and Anderson 1975).Average τ c is proportional to the fault length divided by the rupture velocity and also related to M o (Kanamori and Brodsky 2004).For Japanese earthquakes, Iio (1986) found a change in scaling of τ c versus M o from small to large events.For the aftershocks of the 28 June 1991 Sierra Madre, California, USA earthquakes, Ma and Kanamori (1994) observed a small exponent value of the τ c -M o relation when M o < 10 21 dyne-cm and a large exponent value (~1/3) when M o > 10 21 dyne-cm.Kanamori and Brodsky (2004) observed M o ~ τ c 3 in a large range of M o .
Using the two-point correlation method, Kagan and Knopoff (1980) determined the distribution of distances between pairs of both epicenters and hypocenters of earthquakes.Their results show that the number of events per unit volume at a distance, R, from any earthquake is proportional to R h -, where η approximates 1 for shallow earthquakes and increases to 1.5 or possibly higher for deeper ones.This indicates an existence of fractal characterization of the spatial distribution of earthquakes.Kagan and Jackson (1991) reported a power-law temporal distribution of long-term variation in seismicity.Wang andLee (1995, 1997) reported that the earthquake sequences show multifractal behavior.A detailed description concerning scaling laws of earthquakes can be found elsewhere (e.g., Scholz 1990;Turcotte 1992;Kanamori 1994).
It is necessary to study the rupture processes of earthquakes to understand the physics of their scaling laws.Such processes are very complicated, however, and cannot be completely solved using a simple model.Several factors must be taken into account for modeling.A minimal set of ingredients includes: plate tectonics, brittle-ductile fracture rheology, the stress re-distribution after fracture, friction, the geometry of faults, and the healing process from the dynamic friction strength to the static friction strength after a fault stops moving.In addition, stress corrosion (Anderson and Grew 1977;Atkinson 1984) is often considered to be a mechanism for fractures.Since these ingredients are only partly understood, a comprehensive equation to describe fault dynamics has not yet been established.Nevertheless, some significant models have been developed to approach fault dynamics.
Burridge and Knopoff (1967) proposed a one-dimensional dynamical lattice model (abbreviated as the 1-D BK model henceforth) to approach fault dynamics.This model and its revised versions have been widely applied to simulate the occurrences of earthquakes (Otsuka 1972;Yamashita 1976;Rundle and Jackson 1977;Cohen 1979;Cao andAki 1984/85 and1986;Carlson and Langer 1989a, b;Carlson 1991a;Carlson et al. 1991;Wang 1991Wang , 1993Wang , 1994Wang , 1995aWang , b, 1996Wang , 1997a, b;, b;Knopoff et al. 1992;Shaw et al. 1992;Shaw 1993Shaw , 1994Shaw , 1995;;Pepke et al. 1994;Xu and Knopoff 1994;Schmittbuhl et al. 1996).Before 1989, the number of mass elements of the model used was small, and, thus, simulation results were not good enough.Since 1989, the number of mass elements of the model has been largely increased for obtaining good results because of an increase in memory and computing speed of computers.Most of the studies were based on the 1-D BK model, and only a few (e.g., Calrson 1991b;Miyatake 1992;Nakanishi 1992;Chen 1996;Wang 2000) were based on the 2-D model.In addition, viscous force was also introduced by some authors (Burridge and Knopoff 1967;Rundle and Jackson 1977;Pelletier 2000;Wang 2007) to study aftershocks and source rupture.
Several authors (Carlson et al. 1994;Wang et al. 1995;Wang 1999;and Pelletier 2000) have individually published review articles concerning studies of seismicity based on the 1-D BK model.The first group of authors focused their attention on studies of the intrinsic properties and FM relation based on the BK model.The second and third authors mainly reviewed studies of the FM relation.The fourth author paid attention to studies of structural heterogeneity and viscous effect.It seems to be significant to have a more comprehensive review article to show studies of scaling laws of earthquakes based on the 1-D BK model.Since there is no characteristic length in the BK model, this model has been widely used to study self-organized criticality (SOC) of dynamical systems (Bak et al. 1987(Bak et al. , 1988)), including earthquakes.In addition, there are some publications for studies of nonlinear behavior of the 1-D BK model.Although the two topics are quite interesting, they are out of the scope of this article and it will not be included.
Hence, an attempt in this article is made to review publications, especially for those done by Taiwan's scientists, on the intrinsic properties of the model and studies of several scaling laws based on the 1-D BK model.Those scaling laws in this review include: Omori's law, the FM relation, the relation of rupture length versus seismic moment, the relation of source duration time versus seismic moment, the frequency distribution of rupture lengths, and the earthquake source power spectra.In addition, the studies of frictional and viscous effects on earthquake rupture are also reviewed.

description of the Model
The 1-D BK model consists of a chain of N mass elements of equal mass, m, and springs with each mass element being linked by two coil springs of strength, K c , with two nearest neighbors and each mass element also being pulled through a leaf spring of strength, K l , by a moving plate with a constant velocity, V p .This system is illustrated schemati-cally in Fig. 1.Initially, all mass elements rest in an equilibrium state, and the spacing between two mass elements is 'l'.The n-th mass element is located at position u n , measured from its initial equilibrium position, along the x-axis, which is in the direction of motion.Furthermore, the n-th mass element is subjected to a frictional force, F n .The equation of motion of the n-th mass element of the system is: In this equation, v ( du d ) t n n = is the sliding velocity of the n-th mass element.Obviously, the spacing, l, is not an explicit parameter in Eq. (1).Yamashita (1976) compared Eq. ( 1) with a finitedifference equation, which is an approximation of a 2-D plain-strain-type wave equation in the neighborhood of a fault surface.His results led to the relations of K c and K l to two Lame's constants (λ and μ) and the ratio of the S-wave velocity, v s , and P-wave velocity, v p , of the material, i.e., [2 ( . δy and δz are the spacing units along and perpendicular to the axis of the model, respectively.When δy = δz, both z y d d and y z d d are 1, thus leading to the direct relations of K c and K l to physical parameters of the materials.The plate velocity V p is ~10 -12 m s -1 . The equation of motion essentially consists of two processes: The first one is the coupling process between the moving plate and a mass element through the leaf spring L. The other one is the generation of "self-stress" defined by Andrews (1978).Self-stress originates from the joint effect of the coil spring K c between two mass elements and the leaf spring K l .The coil spring K c plays a role only in transferring energy from one mass element to another; thus, it does not change the total energy of the system.However, the spring K l plays two roles: One is to provide energy to the system from the driving force caused by the moving plate, i.e., the K l V p t term in Eq. ( 1), and the other is to take energy from the system.This indicates that the spring K l can change the total energy in the system.Therefore, the stiffness ratio ( K) s K c l = is a significant parameter represent-ing the level of conservation of energy in the system.Larger s shows that the coupling between two mass elements is stronger than that between a mass element and the moving plate.This results in a smaller loss of energy through the K l spring, thus indicating a higher level of conservation of energy in the system.Smaller s indicates a lower level of conservation of energy.When K c is much smaller than K l or K c equals zero, Eq. ( 1) becomes of the following form: The system loses the coupling between the mass elements, and, thus, each mass element slides independently.Hence, the system can only generate a small event consisting of either one single mass element or only a few mass elements.In other words, the system cannot self-organize itself to form a large event consisting of a large number of mass elements.The magnitude of each event is proportional to the amount of force drop.Unless the distribution of the force drop is initially designed to be a power-law function, the resultant FM distribution will not necessarily be a power-law function.There might be two extreme cases of occurrence of events: one of them shows a sequence of equal-magnitude small events when the distribution of the breaking strengths is inhomogeneous.The other displays a very large event consisting of all independent small events when the distribution is homogeneous and all isolated mass elements can slide simultaneously.On the other hand, when K l is much smaller than K c or K l equals zero, the effect due to the leaf spring disappears, and the coupling between two mass elements dominates the behavior of the system.For the case not including the frictional action, the system is mainly a conservatively self-organized one.Numerous authors have studied this kind of system, which can easily show SOC, using the cellular automata iteration.Since the fault system is a dynamic one with dissipation, its stiffness ratio, s, must be a non-zero finite value.
For a slip weakening, single-degree-of-freedom springslider model, the stiffness K of the spring is considered to be the main parameter controlling the instability of the model system (cf.Rice 1979 andLi 1987).For such a model, small K (less than a critical value) rather than large K can produce an unstable rupture.Based on a simple two-dimensional anti-plane strain softening model, Stuart (1981) considered the ratio K K f s , where K f and K s are the stiffness of the fault zone and that of the elastic surroundings, respectively, to be a significant indicator of earthquake instability.He stated that instability occurs when the ratio reaches unity.Stuart (1986) redefined K K s f as the stiffness ratio to indicate the instability of the system.Obviously, this parameter is just the ratio between the stiffness of the fault zone to that of the elastic surroundings and cannot directly display the coupling between them unless a constitution law is given to describe the correlation between the fault and the elastic surroundings.Although the stiffness K was considered in a series of work by Stuart and his co-authors (cf. Stuart 1988), they did not study in depth the effect of the variation in K on the earthquake rupture.In numerous studies through the cellular automaton iteration, the 1/s value was designated to be zero, i.e., K l = 0, (cf.Takayasu and Matsuzaki 1988;Bak and Tang 1989) or to be a very small value (Brown et al. 1991).In those studies, the effect due to the coupling between a mass element and the moving plate was mostly ignored.In this kind of modeling, the system can self-organize itself very easily and shows SOC.However, in some other studies, the effect of the variation of the stiffness ratio on the scaling of seismicity was included.Based on a simulation for a two-degree-of-freedom earthquake model, Nussbaum and Ruina (1987) addressed the importance of the stiffness ratio (called the coupling ratio in their article) on the slip pattern.The value of the stiffness ratio, s, which best models the earth, has not yet been well determined.Based on theoretical considerations, Carlson and Langer (1989a, b) used large values of ~1000.The values of s used by Wang (1991Wang ( , 1993Wang ( , 1994Wang ( , 1995aWang ( , b, 1996Wang ( , 1997a, b;, b;Wang and Huang 2001) in numerical simulations were less than 120.Knopoff et al. (1992) and Xu and Knopoff (1994) used very small values of 0.625 to 2.632 in simulations.
The boundary conditions at the ends of the model affect the numerical results.Christensen and Olami (1992a) first stated that the scaling exponent depends on the boundary condition.However, their result shows that the difference of the scaling exponents for free and open boundary conditions is reduced as the parameter ( ) s s 4 1 + increases.A periodic boundary condition has been applied at the two end mass elements in solving Eq. (1) for most studies.However, Xu and Knopoff (1994) supposed that systems with periodic boundary conditions will ultimately display runways if the computations are extended to sufficiently long times.

Friction
Friction is a very complicated physical process.From laboratory experiments, Dieterich (1972) first observed the time-dependence of the static frictional strength of rocks in laboratory experiments.Dieterich (1979) and Shimamoto (1986) pointed out the velocity-dependence of the dynamic friction strength.Dieterich (1979) and Ruina (1983) proposed empirical velocity-and state-dependent friction laws.Essentially, the velocity-dependent friction law includes two processes: the velocity-weakening process and the velocity-hardening one.A detailed description of the generalized velocity-and state-dependent friction law can be found in Marone (1998).
For numerical simulations and theoretical analyses, several simple friction laws have been used.Burridge and Knopoff (1967) first considered a velocity-dependent, weakening-hardening friction law.Cao and Aki (1984/85) used a displacement hardening-softening friction law.By comparing a slip-dependent friction law and a velocity-dependent one through numerical computations, Cao and Aki (1986) concluded that the two kinds of friction laws cause different effects on simulation results.Carlson and Langer (1989a, b) considered a velocity-weakening friction law: where v 1 is a particular speed that characterizes the velocity dependence of F. The function F(v) ranges between F o at v = 0 and decreases monotonically to zero as v becomes large.Carlson et al. (1991) made two changes from the definition of F(v) as described in Eq. ( 2).First, the breaking strength (i.e., the static friction strength) at v = 0 can have any value in the range (-3 , 1].In other word, F(v) is a multivalued function at v = 0.The frictional force is a negative infinity when v < 0, thus indicating that no backward motion is allowed.The second change is that while the breaking strength remains at unity, slipping friction begins at φ = 1 -σ, where σ is a small value used by Carlson et al. (1991) for the initiation of an event.This quantity can affect the lower bound magnitude of localized events as mentioned below.Hence, the revised frictional law used by Carlson et al. (1991) has the following form: In the normalized, non-dimensional form of the equation of motion used by Carlson and her co-authors, F o in Eqs.
For the first-order approximation, Wang and Knopoff (1991) considered a piece-wise, linearly velocity-dependent weakening-hardening friction law (as shown in Fig. 2) in the form: which is a simplified form of the friction law proposed by Burridge and Knopoff (1967).F o denotes the breaking strength.The decreasing rate, r, and increasing rate, γ, of dynamic friction strength with sliding velocity are two parameters of the model.As shown in Fig. 2, Eq. ( 4) is defined only for v 0 F and F(v) is a negative infinity when v < 0. This means that no backward motion is allowed.When v = v c , F(v) is the minimum value, gF o .Wang (1995a) called g the frictional drop ratio.The value of g is positive yet smaller than 1.Smaller g produces a larger force drop, thus being able to generate a larger event.Hence, a drop in the frictional strength from F o to gF o behaves like a source supplying additional energy to a mass element for sliding.From analytic analysis, Wang (2002) deeply studied the commonly used velocity-and state-dependent friction law and claimed that Eq. ( 4) is a good approximation of it.
The distribution of the breaking strengths affects dynamical behavior of the system.Carlson and her co-workers used an almost uniform distribution of the breaking strengths.Thus, de-localized events, for which all mass elements of the model are in an unstable state, can be easily generated in their simulations.Nussbaum and Ruina (1987) claimed that such a homogeneous fault stress is generally unstable.It is known that the fault zones where earthquakes occur are usually quite complicated.Seismological and geological observations show that the mechanical properties and geometry of a fault zone are heterogeneous (cf.Wang 2006a, b).From laboratory experiments, Mogi (1963) addressed the importance of the heterogeneity of the materials of the fault plane on the variations for seismicity and the b-value.But, in contrast, based on laboratory results, Scholz (1968) stressed that the state of stress, rather than the heterogeneity of the material, plays the most important role in determining the b-value.The breaking strength is the main mechanical parameter representing the state of stress over the fault zone for generating ruptures and is one of the most important properties certain to influence seismicity and its scaling laws.Das and Aki (1977) and Aki (1979) defined a barrier model and Kanamori and Stewart (1978) defined an asperity model to describe such an inhomogeneous distribution of the breaking strengths of the fault zone for earthquake occurrences.Based on a single rider model, Nur (1978) used a one-body dynamical model to study the effect of displacement-dependent or position-dependent friction on a rupture.His results reveal the importance of the inhomogeneities of the breaking strengths over the fault plane on the propagation of a rupture.He related the rupture velocity to the gradient of the breaking strengths.Rice (1993) and Knopoff et al. (1992) also stressed the importance of inhomogeneous friction on earthquake ruptures.In addition, a reason to form de-localized events might be the use of uniform friction for which the fault plane can break in a short time.
In earlier studies, several authors (e.g., Yamashita 1976;Rundle and Jackson 1977;andCao andAki 1984/85, 1986) used an inhomogeneous distribution of the breaking strengths in seismicity simulations.In practice, different distribution functions can be used to describe the inhomogeneities.Field survey results (Scholz and Aviles 1986;Aviles et al. 1987;Okubo and Aki 1987) and laboratory observations (Brown and Scholz 1985) have suggested that the geophysical and geometrical properties over the fault planes have, in general, a fractal distribution.Fractal properties are commonly found in natural phenomena (cf.Mandelbrot 1982;and Turcotte 1992).Mandelbrot (1982) defined the fractal dimension, D, to describe the fractal set of objects or a fractal geometrical structure.Wang and Knopoff (1991) first suggested a fractal distribution of the breaking strengths for numerical simulations.Since a fractal distribution of a physical property leads to nonlinear behavior, the use of such a fractal distribution makes the model become nonlinear.Wang (1996) applied an infinite BK model together with a linear velocity-weakening frictional law to deduce three types of the propagation of motion of mass elements.Such three types are related to three kinds of velocityweakening friction, which depend on three parameters of the friction law, i.e., the decreasing rate r, the stiffness K l , and the mass m.The three kinds of friction are (1) subsonic-type friction when r > 2 (K l m) 1/2 , (2) sonic-type friction when r = 2 (K l m) 1/2 , and (3) supersonic-type friction when r < 2 (K l m) 1/2 .Obviously, supersonic-type friction leads to non-causal ruptures, because the propagation velocity is greater than the sound speed.Knopoff et al. (1992) stated that when r = 2 (K l m) 1/2 (r was denoted by α in their article.),the system is asymptotic to dispersive-free elasticity in the continuum limit.Wang (1996) also mentioned that the dynamical friction strength with large r or r = 3 behaves like an impulse, and it can supply energy to a mass element in a very short span, thus resulting in a large initial acceleration to the mass element.But, the dynamical friction strength with small r cannot work in this way.Therefore, large r is more capable of generating large events than small r.However, velocity-weakening friction with large r intrinsically prohibits the formation of very large events, for which a large number of mass elements slide in a short span.Such an event was called the 'de-localized event' by Carlson and Langer (1989a, b) and can be very easily generated from the models used by Carlson and her co-authors.When the friction law has the form of Eq. ( 2), the related decreasing , whose value changes from v 1 1 to 0 when v varies from 0 to 3 .Whereas, in Carlson et al. (1991), Eq. ( 3) was used, and, thus, the decreasing rate is , whose value changes from 1 to 0 when v varies from 0 to 3 .Their model basically exhibits supersonic behavior with r < 2 (K l m) 1/2 , for which the abovementioned intrinsic effect does not exist.Hence, their model is potentially capable of producing very large events.
When a mass element stops motion after sliding, the healing process of friction from dynamic friction strength to breaking strength, including the type and delay time of the healing process, would influence the next rupture.In other words, the consequence of non-instantaneous healing must be significant for seismicity (Rundle and Jackson 1977).Cao and Aki (1986) stated that the non-instantaneous healing lengthens the time needed for a fault slip to stop, reduces the interaction between different fault segments and, finally, counteracts the smoothing effect.Wang (1997a) studied the effect on the scaling of seismicity due to the healing process of friction.He defined a parameter to be the ratio of the frictional healing rate, h, to the tectonic loading rate, K l V. From simulation results, he claimed that the ratio h K V l p is only a minor factor in affecting the scaling of seismicity.
Other than these, the effects caused by a non-instantaneous healing process were not included in other papers.

Intrinsic Properties of the Model
Rice (1993) and Rice and Ben-Zion (1996) discussed whether the self-organizing process is capable of generating slip complexity on a spatially uniform fault.They argued that all results from the self-organizing models are an artifact of the use of simplified constitutive laws and also sensitive to the spatial cell size, h, used for numerical simulations.They also suggested a characteristic length-scale h* based on the rate-and state-dependent frictional law proposed by Ruina (1983).They called those models with h > h* the inherently discrete models, for which some GR-type range of small events seems to be generic.For those self-organizing models, they suggested that h* = 0. Xu and Knopoff (1994) stated that slip complexity should not be regarded as a generic feature of the nonlinear dynamics of a smooth fault, whether we model this fault from the framework of continuum mechanics or from dynamical lattice system, but rather should be considered as part of an evolutionary transient process.
They stressed that inhomogeneity is a critical property of earthquake models.On the other hand, Shaw (1994) presented a counter-example to this suggestion.His numerical results are shown to be independent of the spatial discretization for small discretizations compared to the characteristic length-scale.The qualitative features of the complexity produced are seen to be invariant with respect to two very different types of small-scale cutoffs, implying a universality of the results with respect to the details of the small-scale cutoff.From studies of the rupture propagation along the anti-plane fault models with single and twin asperities in the presence of nonlinear velocity-weakening friction, Cochard and Madariaga (1994) observed a complex distribution of stresses after a rupture.Ben-Zion and Rice (1995) and Rice and Ben-Zion (1996) argued that the possible existence of weaknesses in the models used by Shaw (1994) and Cochard and Madariaga (1994).Wang and Huang (2001) considered that the two factors are both important in controlling slip complexity of earthquake faults.They also stressed that only under inhomogeneous breaking strengths, can nonlinear friction play a significant role on slip complexity.

definitions of Model Events and Magnitude scales
For a certain mass element, when the sum of the driving forces due to the moving plate and spring forces from its neighbors exceeds the breaking strength, it is accelerated and starts to slide.After a while, the increase in either the spring force due to the change in the relative positions of the mass element and its neighbors or in the dynamic friction strength with sliding velocity decelerates the motion.Finally, the mass element stops moving and sticks, and the results in a drop in the total force.However, the moving plate, which always loads the mass element, increases the total force on the mass element to reach the breaking strength, and then to push it to slide again.
The displacement of a mass element is measured from its new equilibrium position to the one where it sticks after sliding.This position becomes a new equilibrium position for the next stage of motion.Since several mass elements might slide almost simultaneously within a certain time span, the sum of the displacements of the related mass elements in such a time span provides the time history of the displacement.Such a time history is considered an event.An example to show the space-time patterns (abbreviated as the ST pattern) of synthesized events for four values of s, i.e., 5, 40, 80, and 120, are shown in Fig. 3.The line segment linking up the slid mass elements represents an event.The longer line segment consists of a larger number of mass elements and represents a bigger event.In all cases considered, the number of mass elements slid during an event is about one or two.It is evident that different values of s produce different ST patterns.For small s (for instance 5 and 40), the number of mass elements slid during an event is generally small; in contrast, for large s (for instance 80 and 120), a lot of larger events with longer line segments appear.This indicates that for larger s, a larger number of mass elements can be driven to an unstable state almost simultaneously during a time span, thus leading to a bigger event.In the four cases, the number of events in the range with strong breaking strengths is evidently smaller than that in the range with weak breaking strengths.When s = 40 and 80, the events repeat themselves very frequently at two small ranges with very low breaking strengths.
Since the seismic energy E s is proportional to the maximum slip, u max , and the related force drop Δf, i.e., E f u max s , Wang and Knopoff (1991) defined the logarithmic value of the sum of the seismic energy of mass elements slid during an event to be the magnitude, M, i.e., This magnitude is an energy-based magnitude rather than the commonly used magnitude based on the peak amplitude of the seismogram.This magnitude was exclusively used in Wang's studies.The magnitude used by Carlson and her co-workers was based on the earthquake moment, which is the sum of the displacements of mass elements slid during an event, i.e.,

omori's law for Aftershocks
Only few studies about aftershock activity are based on the 1-D BK model.Burridge and Knopoff (1967) stated that by the introduction of viscosity into the model, aftershocks occur following a major shock, and the time sequence of the number of simulated aftershocks can be fitted by Omori's law.Rundle and Jackson (1977) mentioned that foreshocks and aftershocks were observed as the amount of anelasticity put into the system.Pelletier (2000) constructed a structurally-heterogeneous 1-D dynamical lattice model coupled to a viscous asthenosphere.His model can generate foreshocks and aftershocks which follow Omori's law.

relation of Frequency Versus Magnitude as well as Energy
Gutenburg and Richter (1944) first proposed a frequency-magnitude (FM) law, i.e., log N a bM = -, to describe the earthquakes in southern California.This law appears to hold not only for mainshocks but also for aftershocks and not only for global earthquakes but also for regional ones (cf.Utsu 1961).It also holds for rock fractures (cf.Mogi 1967).The b-value, usually ranging from 0.8 to 1.2, varies from region to region and is also different for various time periods (cf.Wang 1988).In principle, this relation has been taken as an indication of self-similarity or scale-invariance of earthquakes at all magnitudes.However, observations usually show that self-similarity exists only in a range of magnitudes.For individual matured faults, however, the number of large earthquakes and the maximum earthquake usually cannot be interpreted by observed FM relation obtained from small earthquakes (cf.Stirling et al. 1996).The Gutenberg-Richter FM relation implies the number rate of occurrence of earthquakes with energies E in the form: , where p is almost 2/3.If the rate of energy released in earthquakes [ ( ) ] E dN E dt dE # is to be finite, p cannot be greater than 2/3 (cf.Knopoff et al. 1992;Knopoff 1997).Hence, the scale-invariance implied by the GR-type FM relation for the low-energy branch cannot extend to the high-energy one, and there must be a cutoff to the size spectrum at the high-energy end.Knopoff and his co-workers stressed that the low-energy branch is indeed due to a universal process, which does not exhibit SOC.They also argued that the high-energy branch is a response to physical processes that are strongly influenced by the geometry of earthquake fault.The geometry of fault maps varies from region to region.Hence, the high-energy branch must have a local character rather than a universal one, unless the major faults on which large earthquakes occur have some universality of their own geometrical structures.However, this does not seem likely.Pacheco et al. (1992) pointed out the existence of a break in the FM distribution, from small to large earthquakes.The magnitude related to such a break is about 7.5 for shallow earthquakes, about 6.4 for intermediate and deep earthquake, and about 5.9 for the earthquakes on the transform faults in the oceans.Okal and Romanowicz (1994) also reported a similar result, even though their break magnitudes are different from those obtained by Pacheco et al. (1992).It has been hypothesized that the break might occur at a point where the dimension of the event equals the down-dip width of the seismogenic layer.Such a change of the scaling FM relation is comparable with the observations (e.g., Shimazaki 1986) that the seismic moment released in small earthquakes scales differently with rupture length than it does for large events, where the crossover between small and large events is associated with a rupture dimension equal to the down-dip width of the seismogenic layer.Knopoff (1996a) found that there is an absence of such a break for California earthquakes occurring during the 1933 -1992 period.In addition, no such a break can be found from simulation results as described below.From a large set of earthquake source parameters for more than 200 global earthquakes, Wang and Ou (1998) did not find any change in the scaling relation of seismic moment versus fault length from small to large earthquakes.
Early studies of the physical factors in affecting the bvalue were based on laboratory work on rock fractures.Mogi (1967) considered the effect of the degree of heterogeneity of the media on the b-value.Scholz (1968) correlated an increase in the b-value with a decrease in the ambient stress level.Burridge and Knopoff (1967) first studied the FM relation based on their 1-D model.They stated that the scaling exponent of the energy-frequency relation is about 1.According to the simulation results using the 1-D BK model, Rundle and Jackson (1977) found that linear behavior of the GR-type FM relation is not immutable but rather is dependent on the mechanical properties of the faults.From the fragmentation of materials and the fractal distribution of the strain and stress of crustal deformations, Turcotte (1986a, b) interpreted the FM relation.King (1983) considered a geometrical origin of the b-value based on self-similar fault systems.The cellular automatum model, percolation theory, and some other statistical physics models have also been applied to study the FM relation (cf.Chelidze 1986;Sornette et al. 1991;Ito 1992;Lomnitz-Adler 1993;and Main 1996).Based on a quasi-dynamic elastic model, Ben-Zion and Rice (1995) stated that FM statistics are approximately self-similar for small events, with b = 1.2, but strongly enhanced with respect to self-similarity for events larger than a critical size.
From analytic studies and numerical simulations based on a modified version of the 1-D BK model, Carlson and Langer (1989a, b) classified the events into three types: microscopic, localized, and de-localized events (as shown in Fig. 4).Theoretically, they defined several parameters, i.e., K m Only localized events in the magnitude range of from M 1 to M 2 exhibit a GR-type FM relation.Since the real fault zones are, at least, twodimensional, their analytical results could only partly interpret observations.The number of microscopic events in the magnitude range from M s to M 1 is remarkably increased with magnitude, and they do not follow a power-law function.In the regime of de-localized events in the magnitude range from M 2 to M d , there is a pronounced peak in the FM distribution, and, thus, it is impossible to extrapolate the magnitudes of de-localized events from the FM relation of localized events.The general feature of simulated FM distribution obtained by Carlson and her co-authors is shown in Fig. 4.
Based on the 1-D BK model with inhomogeneous breaking strengths, Knopoff et al. (1992) found that the near-saturation of small events is identified with finite lattice spacing effects and the roll-off of large events is associated with the constraint that all fractures are confined, and hence that there must be a maximum event.They also stressed that a self-organizing, spatially localized sequence of events constrained by spatial fluctuations and the usual GR-type FM relation is a correlation of the geometry of localization.In addition, they also found that the linear and the roll-off intervals should be fitted by a gamma function, which span both, because they are caused by the same mechanism.Schmittbuhl et al. (1996) reported that two distinct regimes in the statistical distribution of event sizes and magnitudes are separated by a characteristic size, L*, which depends on the elastic stiffness and the dissipation ratio.A characteristic length, L c , related to L*, controls a crossover between two different dynamical regimes.For events of size smaller than L c , the system exhibits scaling laws.
From the above-mentioned expressions for M 1 and M 2 , the magnitude range of localized events is Carlson et al. 1991).It is obvious that the range of localized events is essentially influenced, at least, by two factors.The first one is the stiffness ratio s (denoted by l 2 in their articles) actually reflects the dissipation of energy through the coupling between the fault and the plate.The second one is the amount of the initial drop in the frictional strength when a mass element begins to slide.Furthermore, ~( ) log M s / 10 3 2 d , because the value of σ is almost a constant and only yields a minor effect on the estimate of δM.From observations, δM is usually in the range of from 2 to 5, thus leading to s = 10 to 10 3 , which were usually used by Carlson and her co-authors in numerical simulations.Aki (1981) Hirata (1989) and Wang and Lee (1996) reported a negative correlation between the two parameters for earthquakes in Japan and Taiwan, respectively.Hirata (1989) stated that Aki's fractal dimension of the geometry of fault planes is a special case of the capacity dimension of asperity or barrier distribution, in which all asperities or barriers are connected to each other without isolation.And in such a distribution, the dimension can be regarded as the fractal dimension of the surface of the fault plane.Yet, this is not necessarily true for the observed seismicity produced from various fault planes.Wang (1991) studied the correlation between the b-value and the fractal dimension of the distribution of breaking strengths from synthetic seismicity based on the 1-D BK model.From the numerical results shown in Fig. 5, he concluded that the b-value is not noticeably dependent upon fractal dimension.This result is different from the theoretical speculation by Aki (1981) and by Turcotte (1986a, b) and from observations (Hirata 1989;Wang and Lee 1996;Chen et al. 2006).A possible reason for the difference might be that Wang (1991) considered a fractal distribution of the breaking strengths of the fault, while the others applied a fractal geometrical structure of the fault plane.
Friction would be a factor in affecting the FM relation.Carlson and Langer (1989b) reported that the b-value is affected by the parameter α as defined above.The b-value first increases with α and then becomes a constant when α is larger than a certain value.Wang (1996) stressed that seismicity patterns and the b-values are different for the three types of friction as mentioned above: The largest b-value is associated with supersonic friction, the intermediate one with sonic friction and the smallest one with subsonic friction.He also stated that large r is more capable of yielding localized events and of prohibiting the generation of de-localized events than small r.The magnitude range of localized events slightly increases with r and the related b-value decreases with increasing r.In simulation results obtained by several authors (Knopoff et al. 1992;Xu and Knopoff 1994;andWang 1994, 1995a), de-localized events, obtained by Carlson and her co-authors (as shown in Fig. 4), cannot be found, even though there are few large events obtained by Wang when r < 3 (i.e., supersonic and sonic friction).The number of microscopic events obtained by Knopoff and co-authors and Wang is smaller than that obtained by Carlson and her co-authors.In the case with r = 3 , the friction immediately drops from breaking strength to dynamic friction strength after a mass element starts to slide.The drop behaves like an impulse to provide additional energy to a mass element for advanced sliding.Immediately after the frictional force drops, the dynamic friction strength again increases with sliding velocity.This decelerates the motion.On the other hand, in the case with r = 1, the work caused by a decrease in the frictional force rises very slowly.Thus, the restoring force caused by a change of the relative distance between the slid mass element and its nearest neighbor increases gradually, thus leading to a resistance to the motion of the mass element.Therefore, the case with r = 1 has less potential than the case with r = 3 to produce large events.Furthermore, Wang (1996) stated that the FM distributions for two different values of g are almost similar for localized events, but different for large events.Smaller g can lead to a greater number of large events than larger g.This is due to the fact that the drop of the friction strength from F o to gF o behaves like an impulse to push an object to slide.Smaller g will provide more energy, thus being capable of resulting in a greater number of large events, than larger g.
Although the magnitude range of localized events obtained from analytic studies by Carlson and her co-authors is proportional to log(s 3/2 ), Carlson et al. (1991) reported that the b-values for three large values of s, i.e., 36, 64, and 144, are almost around 1.0 and independent of s.For the 2-D BK model, Huang et al. (1992) also obtained a scaling exponent of about 1.36 for five values of s, i.e., 10, 15, 20, 30, and 40, using the cellular automaton iteration.On the other hand, simulation results obtained by Xu and Knopoff (1994) and Wang (1994Wang ( , 1995a) ) showed that the b-value decreases with increasing s.The simulation results by Wang (1994Wang ( , 1995a) ) will be discussed in depth below.In addition, Nakanishi (1990) calculated the FM distributions for four small values of s, i.e., 2.00, 2.83, 4.50, and 9.50, using the cellular automaton iteration.His results showed a dependence of b on s, even though he did not calculate the scaling exponent.From the results for the isotropic 2-D models based on the cellular-automaton iteration, Christensen and Olami (1992a, b) stated that the b-value decreases continuously as a function of ( ) s s 4 1 + .The b-values are almost constant for large s, but it does decrease with increasing s for small s.
The plots for six values of s, i.e., 10, 30, 50, 70, 90, and 110, are shown in Fig. 5 (Wang 1995a).The number of events used for each case is greater than 6000.The magnitude range, within which the GR-type FM relation exists, seems to decreases with decreasing s.And, it seems more appropriate to interpret the data points by using two regression lines or a curve.It can be seen that the point with M = 0 (denoted as M c ) is almost the point, where all the distributions of logN versus M for different values of s intersect one another.When M < M c , logN decreases somewhat with increasing s and the difference between two values of logN for two values of s at a certain M is small.When M > M c , logN increases with s.Nevertheless, the distribution of logN versus M for the cases with 30 100 s E E are somewhat close to each other.Results indicate that a change of s causes the opposite effects on rupture for large and small events.Smaller s represents a weaker coupling between two mass elements and can only make a smaller number of mass elements slide, thus leading to a larger number of smaller events and a smaller number of bigger events.Whereas, larger s indicates a stronger coupling between mass elements and can induce a larger number of mass elements to slide almost simultaneously, thus resulting in a larger number of bigger events and a smaller number of smaller events.Results show that the lower bound of the magnitude range is almost constant, and both the upper bound magnitude of localized events and the maximum magnitude somewhat increases with s.Simulation results seem to be consistent with the theoretical result obtained by Carlson and Langer (1989a, b).Wang (1994Wang ( , 1995a) ) explored the correlation between the b-value and the stiffness ratio s.The plots of b (based on the cumulative FM relation) versus s are shown in Fig. 6.Wang (1995a) stated that the b-value of the cumulative FM relation is in general smaller than that of the discrete FM relation.This is similar to the conclusion obtained by Main (1992) from observations.The data points for s > 20 are distributed around a line very well, while those for s < 20 depart from the linear trend.The slope value of the regression line for the data points for s > 20 is about -2/3.In Fig. 6, the related regression line is denoted with a solid line.A similar result can also be obtained for the b-s relation based on the discrete FM relation.For the discrete FM distribution, the related slope value is about -1/2.Obviously, there is a power-law function between b and s: b ~ s -2/3 for the cumulative frequency and b ~ s -1/2 for the discrete frequency.Included also in Fig. 6 are the data points for r = 1 (shown by triangles) and those for g = 0.6 (denoted by diamonds) when s = 50 and 100.Although those data points are distributed above the distributions of data points for the cases with r " 3 , they are somewhat distributed along a line being parallel with the solid one.This phenomenon also exists for the discrete frequency distribution.Consequently, for different model parameters b still relates to s in the form of b ~ s -2/3 for the cumulative frequency and of b ~ s -1/2 for the discrete frequency.This seems to suggest that for a certain value of s, the b-value depends upon r and g, but the scaling relations of b versus s are quite robust.Wang's scaling correlation between b versus s is consistent with that shown in Xu and  Knopoff (1994), but different from that shown in Carlson and Langer (1989a, b) and Carlson et al. (1991).The difference between the results obtained by the former authors and those obtained by the latter authors might be due to the use of different distributions of the breaking strengths and friction laws by the two groups of authors.

relation of source rupture duration Versus Earthquake Moment
The earthquake source rupture duration τ c is the time span of the rupture process and is of the order of the dimension of the rupture zone, ξ R , divided by the shear-wave velocity, β, (Brune 1970).From Haskell's model, Brune (1970Brune ( , 1971) related τ c to the corner frequency f c , in the form: f 1 c c x , in the displacement spectrum S(f), where f is the frequency.When f < f c , ( ) S f is approximately independent of frequency, but when f > f c , ( ) S f bends over and decreases as a power function of f, roughly ( ) S f ~ f --2 .Sato and Hirasawa (1973) pointed out that for the circular fault, the width of the displacement P pulse (denoted by t p ) is almost equal to the τ c calculated by v 2 R r p , where the v r is the rupture velocity.Boatwright (1980) approximated t p to τ c , which is calculated by v 1 r because of the consideration of the healing stage, by multiplying a factor.
For Japanese earthquakes, Iio (1986) correlated the radiated energy with the period of the first cycle for P waves (E sP ) and S waves (E sS ) in a power-law form for t p and t s , which is the width of the displacement S pulse, of from 10 -4 to 10 -1 seconds.It is necessary to use two power-law equations to describe two parts of the log-log plots of E sP versus t p and those of E sS versus t s .The first part for the P waves with t p < 1 sec as well as that for the S waves with t s < 10 -1/2 sec is described by a regression line with a scaling exponent value of 4. The second one for the P waves with t p > 1 sec as well as that for the S waves with t s > 10 -1/2 sec is described by a regression line with a scaling exponent value of 3. In other words, two different relations exist for the E sP -t p and E sS -t s scaling and change at the 1-second period of the first cycle or at the radiated energy of about 10 15 erg, the related seismic moment of which is about 10 -21 dyne-cm.
From the broadband seismograms recorded on the TERRAscope for the aftershocks of the 28 June 1991 Sierra Madre, California, USA earthquakes, Ma and Kanamori (1994) calculated the values of τ c and M o .Their results show that it is impossible to describe the log-log plot of τ c versus M o by using a single linear equation.They tried the exponential-law functions with different values of stress drop and attenuation factor to describe the plot.Of course, their selection is one of the ways to describe the data points.An alternative selection is the use of two power-law functions: one with a small exponent value for the part with M o < 10 21 dyne-cm and the other with a large exponent value (~1/3) for the other with M o > 10 21 dyne-cm.Kanamori and Brodsky (2004) observed M o ~ τ c 3 in a large range of M o .The seismic moment, at which the τ c -M o scaling changes, is almost the same for Japanese earthquakes and the 1991 Sierra Madre earthquake sequence, even though their τ c -M o scaling laws in the two regimes of small and large events are different.
The change in the τ c -M o scaling from the regime of small events with a smaller scaling exponent to the regime of large events with a large scaling exponent is important for numerous seismological studies.For instance, the corner frequency f c as well as the related τ c must be given for the prediction of the strong-ground motions for an impeding large earthquake.In generally, the two quantities can be estimated only from the τ c -M o relation obtained from small earthquakes.However, the value of τ c estimated from small events might not be appropriate for large earthquakes.Consequently, it is necessary to study this problem more profoundly.Carlson and Langer (1989a, b) defined an earthquake moment, M', to be the total displacement of a connected set of mass elements, which slide during an event.This definition is different from the original definition of seismic moment, i.e., M o = μδA, where μ, δ, and A are the rigidity, the average displacement, and the source area, respectively.The definition of M' = / u i leads to M' = Nδ, where N is the number of mass elements slid during an event.This leads to M' ~ L, where L is the rupture length of an event and is almost equal to Nl, where "l" is the spacing between two mass elements in the equilibrium state as mentioned above.Since A = LW, where W is the fault width, we have M' ~ δLW ~ δA, because the width, W, of the 1-D model fault can be considered to be unity.This indicates that there is a positive correlation between M o and M'.From numerical results (see Fig. 7 as an example), Wang (1993) observed that the stiffness ratio, s, is a major parameter, while the decreasing rate, r, fractal dimension, D, and roughness, R, are three minor parameters affecting the τ c -M' relation.He also reported that when s > 60, the scaling relations are τ c ~ M' 1/3 for small events and τ c ~ M' 3/5 for large ones.On the other hand, when s < 60, the scaling relation does not change too much for small and large events, and having a form of τ c ~ M' 1/3 .This simulated τ c -M o relation is consistent with the observed one, especially for larger-sized earthquakes.Considering a circular crack model, the seismic moment is given by M o = μπR 2 δ, where μ is the shear modulus, δ is the average fault slip, and R is the radius of the fault plane.The static stress drop during an earthquake is proportional to the displacement divided by a fault dimension, i.e., ~R M R . The source duration time τ c is proportional to R v s assuming that the rupture velocity is proportional to the shear velocity v s .Since . Hence, the source duration time scales with seismic moment in a form of τ c ~ M o 1/3 , if v s and Δσ are almost constant.This indicates that Wang's results for s < 60 are similar to those for a circular crack model.

relation of rupture length Versus Earthquake Moment
A positive correlation between seismic moment (M o ) and rupture length (L) of the faults has been studied for a long time.According to the crack model, Kanamori and Anderson (1975) reported a power-law function in a form of M o ~ L 3 to correlate the two parameters.This M o -L correlation is commonly considered to hold for a wide range of events (see Hanks 1977), with L varying from meters to hundreds of kilometers.Actually, there is a large scatter in the data points of M o versus L between the stress drops from 1 bar to 100 bars, and there may be systematic deviations in the case of very small earthquakes.A close examination of Hank's plots of M o versus L seems to show that the scaling laws are different in two magnitude ranges: one with a larger scaling exponent value (> 3) for smaller events and the other with a smaller one (about 3) for larger events.For Japanese earthquakes, Shimazaki (1986) found that a change from M o ~ L 3 to M o ~ L 2 scaling occurs almost at a point with M o = 7.5 × 10 25 dyne-cm and a fault length of 17 km, which is nearly the thickness of seismogenic layer in the region.Scholz et al. (1986) stated that large intraplate earthquakes consistently have greater moments per unit fault length than interplate events, the difference being about a factor of five.The slip rate of the fault is larger for interplate earthquakes than for intraplate ones.Aki (1992) showed a deviation of the M o -L correlation for the California earthquakes from the average relation in Japan.However, the deviation is larger for strikeslip earthquakes and smaller for thrust ones.Romanowicz (1992) stressed different scaling for strike-slip earthquakes and pure thrust and normal events.She also reported two different scaling laws for small and large strike-slip events, separating at M o = 10 27 dyne-cm: M o ~ L 1/2 for small events and M o ~ L for large ones.However, from more than 200 global earthquakes, Wang and Ou (1998) did not find such a change of scaling, and they also reported that M o ~ L 2 in a large range of M o .Carlson et al. (1991) first studied the correlation between the two parameters using a revised version of the 1-D BK model with a velocity-weakening friction law as described in Eq. ( 2) and with an almost uniform distribution of the breaking strengths.Their results showed that the simulated M' -L distribution cannot be fitted by a simple power law but varies with the size of events: M' ~ L 3/2 for small events, and M' ~ L for large ones.For intermediatesize events, it is not possible to describe the simulated M' ~ L distribution by using a simple and/or single power-law function.As a result, the M' is not given by any single power law of L throughout the entire range of fault lengths.Based on the 1-D BK model, Wang (1995b) studied the M' -L relation.His definition of seismic moment is similar to that used by Carlson et al. (1991).His results show that the process of velocity-weakening friction from breaking strength to dynamic one obviously affects the relation.Only the rapidly weakening process can produce a well-defined powerlaw M' -L relation.The stiffness ratio s is not a significant factor in affecting the scaling relation.Wang (1995b) stressed that the stiffness ratio is a parameter representing the level of energy conservation or the degree of dissipation of energy.Hence, the independence of the M' -L relation on the stiffness ratio shows that the degree of dissipation of energy cannot change this scaling law.Wang's simulated M' -L relations have the following forms: M' ~ L 2 for small events and of M' ~ L for large ones (an example for the simulated M' -L distribution is shown in Fig. 8).However, no transition zone is recognized.Based on continuum models associated with the 1-D BK model, several authors (e.g., Langer et al. 1996;Myers et al. 1996;and Shaw 1997) stated that M' ~ L 2 for the smallest events, and M'~ L for the largest ones.Meanwhile, there is a transition range, where no power-law function can be derived, for the intermediatesized events.The M'~ L relations obtained from simulation results by different authors are obviously different from observed ones.This might be due to a fact that simulations are made based on the 1-D model, while the observations come from natural 2-D fault zones.Hence, it would be significant to study further this problem using a 2-D BK model.

Frequency distribution of rupture length
Field observations show that the fault populations obey fractal geometry in the form of N ~ L -d , where N is the cumulative number of faults having lengths ≥ L and d is the scaling exponent (cf.Davy 1993).A variety of scal-ing exponent values have been observed for fault systems, e.g., 0.89 for the fractures on the Reykjanes, southwest Iceland (Gudmundsson 1987), 2.1 for intraplate faults in Japan (Scholz and Cowie 1990), 1.8 for several data sets obtained by Marrett and Allmendinger (1991), 1.91 for Basin and Range faults and 1.24 at Yucca Mountain, Nevada (Marrett 1994) and 1.3 for the Volcanic Tableland faults (Scholz et al. 1993).For earthquakes occurring in the Geysers geothermal field in northeast California, Sahimi et al. (1993) obtained a value of about 1.9.From laboratory experiments to simulate continental collisions, Sornette et al. (1991Sornette et al. ( , 1993) ) and Davy et al. (1990) reported that the value of d is in the range of from 0.7 to 1.6.However, Davy (1993) and Davy et al. (1995) suggested that a Gamma function in the form of ( ) , where C is a constant, is more appropriate than other functions to fit the discrete frequency-length (FL) distribution obtained from the San Andreas fault system and those from experimental results.They also stated that the characteristic length L o is close to the thickness of the brittle crust.In those observations, the FL power-law holds only in a small range of an order of magnitude in length.
Based on the 1-D BK model, Wang (1997b) simulated the population of earthquake faults.An example of simulated FL distribution is shown in Fig. 9 [L was denoted by ∆ in Wang (1997b)].His results show that the fractal dimension, D, of the distribution and the maximum value of the breaking strengths is a minor parameter affecting the FL distribution, especially for small and intermediate-size events.Nevertheless, larger D can produce longer events than smaller D. The stiffness ratio, s, makes some effects on the FL distribution.The decreasing rate, r, and the increasing rate, γ, and the frictional drop ratio, g, of Eq. ( 4) are three parameters significantly affecting the FL distribution.
Large r, with small γ is more appropriate for generating a power-law FL relation than small r with large γ.Different values of g result in different scaling exponent values of the FL relation existing in an almost same length range.Smaller g (with a larger frictional drop) can lead to a smaller scaling exponent value than larger g (with a smaller frictional drop).When other model parameters are fixed, the scaling exponent value is 1.5 for g = 0.8 and 1.0 for g = 0.6.These simulated scaling exponent values are comparable with the observed ones as mentioned above.

Earthquake source spectra
The body-wave seismic spectrum is controlled by the seismic moment M o and the corner frequency, f c , which is associated with the source dimension.The generally accepted earthquake source functions have either f -2 or f -3 high-frequency spectral decay, and are commonly referred to as ω-square and ω-cubic models (ω = 2πf) (cf.Aki 1967Aki , 1972)).Some authors (cf.Boatwright 1978;Dysart et al. 1988;Patane et al. 1997) claimed that neither of them is appropriate to describe the observations.
Power spectral density with a form of f -2 or f -3 shows a general type of f -β signal.The f -2 signal is considered to be a result of the Brownian motions.Bak et al. (1987Bak et al. ( , 1988) ) proposed self-organized criticality to explain f -β signal.Frankel (1991) assumed that the high frequency energy of a mainshock is produced by a self-similar distribution of subev-  ents, where the number of subevents with radii greater than R is proportional to R -D , D being the fractal dimension.In his model, an earthquake is composed of a hierarchical set of smaller events.The static stress drop is parameterized to R ν , and strength is assumed to be proportional to static stress drop.He found that a distribution of subevents with D = 2 and stress drop independent of seismic moment (ν = 0) produces a mainshock with an f -2 falloff, if the subevents areas fill the rupture area of the mainshock.Based on an ideal system under external random forces, Koyama and Hara (1992) studied the dynamical process of random activation.In their model, the time evolution of the system is described by the Langevin equation, and a scaling rule (represented by an auto-correlation function) to describe the random activation is introduced to generalize the system.Their generalized system predicts the fractional power spectrum f -β (f being the frequency) from a white spectrum to a Lorentzian.Their results show that the exponent β is a function of the fractal dimension of the scaling rule.The fractal dimensions of 2, 1, and of about 0.47 indicate a Lorentz spectrum, an f -1 spectrum, and a power spectrum of f -1.53 type, respectively.
It is significant to study the scaling of earthquake source spectra based on a dynamical model.However, presently only Shaw (1993) has studied this problem.Based on a modified version of the 1-D BK model in the presence of velocity-weakening friction as used by Carlson and Langer (1989a, b), Shaw (1993) obtained the theoretical moment spectra, P(ω).He mentioned that large and small events show different spectra.For large events with M > M*, where M* = 2 a , there are different power-law relations in three angular-frequency regions: In Eq. ( 5), M is the moment, defined as being the sum of the change in the displacements of all mass elements and 2 ( ) ln l 4 2 g v a = . The definitions of α and σ are given previously.When σ is small and α > 1, the exponent ε has almost a value of 2.5.Eq. ( 5) shows that there are two turning points in the moment power spectrum for large events.At low frequencies, the theoretical result is similar to that mentioned by Aki.At medium frequencies, the theoretical power spectrum shows the so-called 1/f noise (cf.Bak et al. 1987Bak et al. , 1988)).At high frequencies, the theoretical result is somewhat between ω -2 and ω -3 models, because ε is about 2.5.This is somewhat different from that shown in Aki (1967Aki ( , 1972)).For large events, Shaw's theoretical source power spectra are more complicated than those obtained from Aki's model.
On the other hand, for small events with M < M*, there are two power-law relations in two angular frequency ranges: In Eq. ( 6), L is the rupture length of an event.Eq. ( 6) exhibits that there is only a turning point, which is dependent upon L, in the moment power spectrum for small events.
Eqs. ( 5) -( 6) show that at low frequencies, M(ω) is almost constant for both small and large events, while at medium and high frequencies, the spectral scaling laws for the two kinds of events are different.Moreover, Shaw's theoretical result for small events is similar to the ω-square model proposed by Aki based on the dislocation theory.Shaw's theoretical results do not seem able to confirm which model, i.e., the ω-square model or the ω-cubic one, is more acceptable than the other for interpreting the source power spectra of earthquakes, especially for large events.The Shaw used friction law was a purely velocity-weakening friction law and an almost uniform distribution of the breaking strengths for numerical simulations.It would be significant to use a velocity-and state-dependent friction law and a heterogeneous distribution of the breaking strengths for further studies.Like the models used by Carlson and her co-authors as mentioned above, Shaw's model was basically a non-causal one.Whereas, Aki (1967) drove the source power spectra based on a causal dislocation model.Therefore, the difference in the models could cause in-consistence between their source power spectra.In addition, Shaw's simulations were based on a one-dimensional model.However, the real earthquake fault must be of two dimensions; and hence, to study this problem, a two-dimensional model is required.

Frictional and Viscous Effects on Earthquake rupture
On 20 September 1999 (at 1747 UTC) an M7.6 earthquake ruptured the Chelungpu fault in central Taiwan (cf.Ma et al. 1999;Shin 2000).The earthquake initiated at 23.853°N, 120.816°E, with a focal depth of ~8 km.The fault is mainly a transpressive one, striking N5°E and dipping 34° to east.Various observations lead to a difference between the northern and southern segments.Wang (2006a) made a detailed review of the differences.Wang (2003) applied a one-body model in the presence of a conventional static/kinetic friction law to approach the motions of each segment.His results suggest that the average displacement on a ruptured area is capable of representing its behavior, which would consist of several asperities with different dimensions, while the predominant period of displacement waveforms is only able to display the oscillations of the major asperity.However, Wang's (2003) simplified model cannot interpret the differences in velocities and accelerations between the northern and southern segments.
In order to understand the detailed properties of faulting, a more comprehensive model is needed.Wang (2007) assumed that viscosity is also a significant factor; and thus, applied a strike-slip-type, two-body BK dynamical model in the presence of both friction and viscosity to approximate the rupture processes of an earthquake along the faultstriking direction.Results show that in addition to friction, viscosity is also an important factor in controlling rupture.For the Chi-Chi earthquake, his simulation results from the model can explain the differences in displacement, velocity, acceleration, and predominant period between the two fault segments.

suMMAry
Studies on the effects on several scaling laws of earthquakes due to model parameters based on the one-dimensional Burridge-Knopoff dynamical lattice model (Burridge and Knopoff 1967) have been reviewed.The scaling laws in utilized in the reviewed studies include: Omori's law, the Gutenberg-Richter-type magnitude-frequency (or energyfrequency relation), the relation between the source duration time and seismic moment, the relation between the source rupture length and seismic moment, the frequency-length relation, and the source spectra.In addition, the effects on scaling exponent, i.e., the b-value of the Gutenberg-Richter-type frequency-magnitude relation are also taken into account.The frequency-magnitude relation as well as the b-value is obviously affected by the changing rate of dynamic friction strength varying with sliding velocity and the stiffness ratio of the model.However, the fractal dimension of the distribution of the breaking strengths is not a significant parameter affecting the b-value.Except for the source power spectra, simulation results are comparable with observations to some extent.This seems to indicate that it is appropriate to apply the 1-D BK model to study the scaling relations between source parameters, even though the values of some model parameters are not yet well known.Consequently, the one-dimensional BK dynamical lattice model acceptably approaches fault dynamics.Of course, a two-dimensional model must be better than a one-dimensional one.Carlson and Langer (1989a, b) for the triggering of an event

Fig. 2 .
Fig. 2. A piece-wise, linearly velocity-dependent frictional law: F o = the breaking strength; v c = the critical velocity; and g = the frictional drop ratio.
nitude is different from the above-mentioned one.Thus, for the energy-based magnitude, the scaling exponent of logN versus M must be similar to 'B' in the relation of N ~ E -B as mentioned previously and different from 'b' in the GR-type FM relation and from 'b' in the relation used by Carlson and her co-authors.But, for simplicity, the notation 'b'will hereafter be used to express the scaling exponent of the FM relation.
α is the ratio of the largest slipping speed to the characteristic speed.

Fig. 4 .
Fig. 4. Figure shows the general pattern of simulated distribution of frequency versus magnitude by Carlson and her co-authors.
postulated a positive relation between the b-value and the fractal dimension, D, in the form of D is the slope of log moment versus magnitude relation, and c is about 1.5.The theoretical relation between the two parameters proposed by Turcotte (1986a, b) are, respectively, b D 3 = and b D 2 = based on different models.For the aftershocks of the 1999 Chi-Chi, Taiwan, earthquake, Chen et al. (2006) observed a positive correlation between b and D. However,

Fig. 6 .
Fig. 6.The log-log plot of b versus s when r = 3 (denoted by open circles).The slope value of the solid line is -2/3.The circle with a cross inside shows the data point when r = 1 (from Wang 1995a).

Fig. 7 .
Fig. 7.The log-log plots of τ c versus M for the model with D = 1.5, s = 100, r = 3 , and g = 0.8.Two solid lines represent the regression lines with different slope values (from Wang 1993).
of the correlation between the number of events per unit volume at a distance, R, and R, i.e., N ~ R breaking strength used by They derived the upper bound magnitude of localized events, i.e.,