A dynamical study of frictional effect on scaling of earthquake source displacement spectra

The scaling of earthquake source displacement spectra is analytically studied based on the continuous form of one-dimensional dynamical spring-slider model in the presence of either linearly slip-weakening friction or linearly velocity-weakening friction. The main parameters of the model are the natural angular frequency, ~o, and the (dimensionless) decreasing rate, D, of friction with slip (or the characteristic displacement) for slip-weakening friction as well as the (dimensionless) decreasing rate, y, of friction with velocity (or the characteristic velocity) for velocity-weakening friction. The analytic solution includes the complementary and particular parts. The former shows the travelling wave and the latter denotes vibrations at a site. The complementary solution exhibits ~-1 scaling in the whole range of ~ for both friction laws. For the particular solution, slip-weakening friction results in spectral amplitudes only at three values of ~. For velocity-weakening friction with y>0.5, the log-log plot of spectral amplitude versus ~ exhibits almost ~0 scaling when ~ is lower than the corner angular frequency, ~c, which is independent on y and increases with ~o. When ~>~c, the spectral amplitude monotonically decreases with ~ following a line with a slope value of −1, which is the scaling ex-


Introduction
The body-wave seismic spectrum, P(~), where ~=2rf is the angular frequency and f is the frequency, is controlled by the seismic moment M o and the corner angular frequency ~c, which is associated with the source dimension.The generally accepted earthquake source functions have either ~-2 or ~-3 high-frequency spectral decay, and are commonly referred to as ~-square and ~-cubic models [Haskell 1966, Aki 1967, Brune 1970, Aki 1972, Kanamori and Anderson 1975].Of course, some authors [cf.Boatwright 1978, Dysart et al. 1988, Patanè et al. 1997] claimed that neither of them is appropriate for describing the observations.Kim et al. [1997] considered the effect of seismic attenuation on the seismograms from which the source displacement spectrum is retrieved.They also suggested a method to separate the attenuation and source effects.Huang and Wang [2002] studied the scaling law of the displacement spectra from the seismograms recorded at nine near-fault stations generated by the 1999 Chi-Chi, Taiwan, earthquake.Results show that the values of corner frequency f c at the nine stations are almost 0.2 Hz.The pattern of displacement spectrum at each station is similar to the theoretical one proposed by Aki [1967].The spectral amplitude is almost constant (or ~0 scaling) when f <0.2 Hz and decays at a higher decreasing rate when f >3 Hz.In the frequency range of 0.2-3 Hz, P(~)~~-b , b varies from 1.63 to 3.04 at the nine nearfault stations, and decreases from north to south.Such a variation might be mainly due to the source effect.Results seem also to suggest that a ~-square model is appropriate for the single-degree-of-freedom rupture processes on the southern fault plane, and a ~-cubic model for the two-degree-of-freedom ones on the northern fault plane.
Scaling of spectral amplitude with ~is generally specified with a form of ~-2 or ~-3 .This shows a type of f -c 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 -c signal.Frankel [1991] assumed that the high frequency energy of an earthquake is produced by a self-similar distribution of subevents.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 n , 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 (n=0) produces an earthquake with an ~-2 falloff, if the areas of subevents fill the rupture area of the earthquake.Based on an ideal system under external random forces, Koyama and Hara [1993] studied the dynamical process of random activation.They applied the Langevin equation to represent the time evolution of the system and took a scaling rule (represented by an auto-correlation function) to describe the random activation for generalizing the system.Their model predicts the fractional power spectrum f -c from a white spectrum to a Lorentzian.The exponent c 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, a f -1 spectrum, and a power spectrum of f -1.53 type, respectively.Herrero and Bernard [1994] and Bernard et al. [1996] proposed the l- square model to approach the theoretical ~-square model.Hisada [2000] proposed a model modified from the l-square model proposed by Bernard et al. [1996] under three assumptions: (1) the spatial wave-number spectrum of the slip distribution falls off as the inverse of the wave-number square; (2) a Kostrov-type slip velocity model is included; and (3) the incoherent rupture time is introduced to model variable rupture velocities.He claimed that his model can produce the ~-square source displacement spectrum.Aki [1967] described the scaling of earthquake source displacement spectra using the dislocation model based on empirical assumptions about some model parameters.Although this approach is good enough, it is still significant to study such scaling based on a dynamical model.Up to date, only Shaw [1993] studied this problem using a modified version of the one-dimensional dynamical spring-slider mode (abbreviated by the 1-D BK model hereafter) proposed by Burridge and Knopoff [1967].Meanwhile, friction controls the rupture processes of an earthquake [e.g., Nur 1978, Knopoff et al. 1992, Rice 1993, Shaw 1993, Wang 1996, 1997, 2000, 2004, 2006a, 2007, 2008, 2009, 2012, 2016].Shaw [1993] considered velocity-weakening friction, as shown below, suggested by Carlson and Langer [1989].He took a homogeneous spatial distribution of static frictional force in his study.His results show a dif-ference in spectra, P(~), between large events and small ones.For large events with M o >M oc , where M o is the seismic moment and M oc is the characteristic seismic moment, there are different power-law relations in three angular-frequency regions: P(~)~~0 for ~< 2r/M o ; P(~)~~-1 for 2r/M o <~<2r/p; and P(~)~~-c for ~>2r/p, where p =2ln(4l 2 /v)/a.The definitions of a and v can be seen in Shaw [1993].When v is small and a >1, the exponent c has almost a value of 2.5.Obviously, there are two turning points in the source spectrum for large events.At low ~, the theoretical result is similar to that proposed by Aki [1967].At medium ~, the theoretical power spectrum shows the so-called 1/f noise [cf.Bak et al. 1987].At high ~, the theoretical result is somewhat between ~-2 and ~-3 models, because c is about 2.5.This is somewhat different from that proposed by Aki [1967].For large events, Shaw's theoretical source spectra are more complicated than Aki's.For small events with M o < M oc , there are two power-law relations in two angular-frequency ranges: P(~)~~0 for ~< 2r/L and P(~)~~-2 for ~> 2r/L, where L is the rupture length of an event.There is only a turning point, which is dependent upon L. Obviously, Shaw's theoretical result for small events is similar to Aki's ~-square model.Shaw [1993] only used a velocity-weakening friction law for numerical simulations.It seems necessary to explore the effect on the source displacement spectrum due to slip-weakening friction.In this work, I will focus on the analytical study of frictional effect on scaling of earthquake source displacement spectra using the continuous form of the 1-D BK model in the presence of two types of friction, i.e., slip-weakening friction and velocity-weakening friction.

One-dimensional spring-slider model
The 1-D BK model (see Figure 1) consists of N sliders of equal mass, m, and springs with one slider being linked by a coil spring of strength, K c , with the other.Each slider is also pulled by a leaf spring of strength, K l , on a moving plate with a constant velocity, V p .At time t=0, all the sliders rest in the individual equilibrium states.The i-th slider (i = 1, …, N) is located at position X i , measured from its initial equilibrium position, along the horizontal axis.Each slider is subjected to a slipand/or velocity-dependent frictional force.The frictional force is denoted by F i (X i ; V i ), where V i =∂X i /∂t, with a static frictional force, F si , at rest.Elastic strain of each slider gradually accumulates due to the moving plate.Once the elastic force at a slider is greater than F si , the static frictional force drops to the dynamic frictional force, F di , and the slider will move subject to F di .The equation of motion is (1) Wang [1995] defined the ratio l =K c /K l to be the stiffness ratio of the system.It is an important parameter representing the level of conservation of energy in the system.Larger (smaller) l shows stronger (weaker) coupling between two sliders than between a slider and the moving plate.This results in a smaller (bigger) loss of energy through the K l spring, thus indicating a higher (lower) level of conservation of energy in the system.Since the fault system is a dynamically coupling one with dissipation, l must be a non-zero finite value.The plate velocity V p is usually small and in the order of ~10 -12 m/s.It is noted that the 1-D BK model can only generate a longitudinal-wave-type rupture [Wang 1996].However, the result obtained from this model is still significant for understanding the P-wave-type source scaling.In crack mechanics, the load is applied remotely on a 3-D extended fault.For the present model, the load is directly linked to each slider.As mentioned below, when the driving force is higher than static frictional force, only the earthquake rupture is taken into account and the driving force and static frictional force will be cancelled out each other and can be ignored.Hence, only a single instability is considered and therefore the problem can be treated easily.The present advantage is that the closed-form analytical solutions can be obtained.

Velocity-dependent friction and slip-dependent friction
The frictional force between two contact planes is classically considered to drop from static one to dynamic one after the two planes move relatively.In fact, friction is a very complicated physical process.From laboratory experiments, Dieterich [1972] first found time-dependent static frictional strength of rocks in lab-oratory experiments.Dieterich [1979] and Shimamoto [1986] observed velocity-dependent dynamic frictional strengths.Dieterich [1979] and Ruina [1983] proposed empirical, velocity-and state-dependent friction laws.In fact, a large debate related to the friction laws governing earthquake ruptures has been made for a long time.Although this problem is important, it is out of the scope of this article and thus will not be described in details.A detailed description of the generalized velocity-and state-dependent friction law and the debates can be found in several articles [e.g., Marone 1998, Wang 2002, Bizzarri and Cocco 2006c, Wang 2009, Bizzari 2011].
For theoretical analyses and numerical simulations of earthquake ruptures, several simple friction laws have been taken into account.Burridge and Knopoff [1967] first considered a velocity-dependent, weakening-hardening friction law.Carlson and Langer [1989] and Carlson et al. [1991] considered a velocity-weakening friction law: F(v)=F o /(1+v/v c ), where F o is the static frictional force and v c is a characteristic speed to specify the variation in F with velocity.F(v) is F o at v=0 and decreases monotonically from F o to zero as |v| becomes large.Several authors [e.g., Carlson and Langer 1989, Carlson 1991, Carlson et al. 1991, Beeler et al. 2008] theoretically modeled earthquakes by using this velocity-weakening friction law.Cochard and Madariaga [1994] and Madariaga and Cochard [1994] assumed that purely velocity-dependent friction models can lead to unphysical phenomena or mathematically ill-posed problems.Hence, those models are very unstable at low values of the fault slip velocity both during the passage of the rupture front and during the possible slip arrest phase.Moreover, Ohnaka [2000] stressed that purely velocity-dependent friction is in contrast with laboratory evidence that is the friction law is not a one-valued function of velocity.Bizzarri [2011a,b] deeply discussed this point.Nevertheless, in order to obtain a closed-form solution the single-valued velocity-dependent friction law is taken into account.The present study can be regarded as a marginal analysis of the effect of velocity-dependent friction on the scaling of source displacement spectra of earthquakes.Analytic results will help us to understand such an effect.
For the first-order approximation, Wang [1995Wang [ , 1996] ] considered a piece-wise, linearly velocity-dependent weakening-hardening friction which is simplified from the friction law proposed by Burridge and Knopoff [1967].The decreasing (weakening) and increasing (hardening) rates of dynamic friction strength with sliding velocity are two main parameters of this friction law.The two rates are denoted, respectively, by r w and r h .The function F(v) is defined only for v ≥ 0 and no backward motion is allowed.When v=v c , F(v) reaches the minimum value, i.e., gF o (0<g<1).The smaller g is, the larger the force drop is.Thus, smaller g is more capable of generating a larger event than larger g.A drop in the frictional force from F o to gF o behaves like an extra source supplying additional energy to a slider for proceeding sliding.The velocity hardening stage is a fast re-strengthening process associated with some healing of slip.Cao and Aki [1984/85] took a displacement softening-hardening friction law: where d is the fault slip and F o , d o , and c are constants.Cao and Aki [1986] applied rate-and state-dependent friction to controlling the motion of a slider.From a comparison between a slip-dependent friction law and a velocity-dependent one through numerical computations, they concluded that the two kinds of friction laws cause different effects on simulation results.
Friction can also be produced from thermal pressurization [cf.Fialko 2004, Bizzari and Cocco 2006a,b, Rice 2006, Wang 2009, Wang 2013].Rice [2006] proposed two end-members models for thermal pressurization: the adiabatic-undrained-deformation (AUD) model and slip-on-a-plane (SOP) model.He also obtained the shear stress-slip functions, x(d), caused by thermal pressurization: where d c =tch/fK and t, c, h, f, and K are, respectively, the density, heat capacity, thickness, friction coefficient, and undrained pressurization factor of the slip zone, for the AUD model; and , where erfc(z) denotes the complementary error function and L * = 4(t c /K) 2 (a t 1/2 +a h 1/2 ) 2 /f 2 V, where a t , a h , and V are, respectively, the thermal diffusivity, hydraulic diffusivity, and slip rate of the fault zone, for the SOP model.The two parameters d c and L * control the shear stress-slip functions of the individual endmember models.Clearly, the two friction laws exhibit exponentially slip-weakening friction for the AUD model and slip-weakening friction for the SOP model.

Static frictional force and breaking strengths
The distribution of the static frictional force (or the breaking strengths) affects dynamical behavior of a fault system.Carlson and her co-workers used an almost uniform distribution of the breaking strengths in their studies.Thus, de-localized events, for which all sliders 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 related to earthquakes are usually quite complicated.Seismological and geological observations show that the mechanical properties and geometry of a fault zone are heterogeneous.For example, the breaking strengths on the Chelungpu fault, along which the 1999 M s 7.6, Chi-Chi, Taiwan, earthquake ruptured, are heterogeneous [see Wang 2006a,b].The breaking strength is the main mechanical parameter representing the state of stress over the fault zone for generating ruptures and thus influences seismicity and its scaling laws.Two phenomenological models have been long accepted to describe the heterogeneous fault strengths: one is the barrier model proposed by Das and Aki [1977] and Aki [1979] and the other the asperity model suggested by Kanamori and Stewart [1978].

Analytical studies and discussion
In this study, I will analytically study the conditions of generating scaling of source displacement spectra of earthquakes utilizing the 1-D BK model in the presence of either velocity-or slip-weakening friction.As mentioned above, the friction law is quite complicated.In order to perform analytic studies, I take the simplified linear friction laws: the velocity and V c = the characteristic velocity) for velocity-weakening friction and F(X)=F o (1−X/X c ) (X=the slip and X c =the characteristic slip distance) for slipweakening friction.The former and the latter are almost the approximations, respectively, of the friction law used by Carlson and Langer [1989] at low velocities and of that used by Cao and Aki [1986] or those proposed by Rice [2006] at small displacements.It is obvious that F o /V c is the weakening rate with velocity and F o /X c is the weakening rate with slip.Hence, in some sense the marginal analyses of generalized friction laws are made in this study.The characteristic velocity and characteristic slip distance are also the slopes of the respective laws, thus representing the respective decreasing rates.A uniform distribution of the breaking strengths is considered in this study.This means that only steady travelling waves are taken into account. Letting and thus X i =x i +V p t.This makes Equation ( 1) become (2) When a slider slips, V p t and V p are, respectively, much smaller than x i and v i , and thus can be ignored.Thus, Equation (2) becomes (3) Obviously, V p is no longer a parameter of the equation of motion. ; .
Dividing Equation (3) by a unit area, i.e., dydz, where y and z denote the coordinates, respectively, along and normal to the 1-D system, leads to (4) Defining t A =m/dydz, k c = K c /dydz, k l = K l /dydz, and f i =F i /dydz, Equation ( 4) becomes (5) To simplify the mathematical manipulation, the continuous form will be normalized on the basis of some normalization factors.Now, the stiffness ratio be- to be the characteristic slip distance of a slider exerted by a force f o through a spring with strength of k l .Larger f o yields longer D o when k l is fixed.Define ~o=(k l /t A ) 1/2 =(K l /m) 1/2 and x = ~ot.The quantity ~o/2r is the predominant frequency of oscillation of a single slider attached to a leaf spring in the absence of friction and thus a characteristic parameter of the model.This gives T o =2r/~o to be the natural period of oscillation.The ratio D o /V p is the loading time for a leaf spring to stretch enough for overcoming the static frictional force, and thus V p /D o ~o is equivalent to the ratio of the slipping time ~o -1 to the loading time.Obviously, D o and ~o are two basic units to scale, respectively, the spatial coordinates, x i , and time, t.
Defining u i =x i /D o and v i =du i /dt and substituting the above-mentioned normalized factors into Equation (5) lead to (6) To further perform the problem easily, like Carlson and Langer [1989] the equivalent differential equation, i.e., the continuous form, is made from Equation (6).The difference terms in the right side of Equation ( 6) are divided by the square of the unit length, a, between two sliders in the equilibrium state and thus Equation ( 6) becomes (7) It is noted that the length 'a' is not an explicit one of the original equation.The limits of la 2 , u i , and v i , are represented, respectively, by h 2 , u, and v as the value of 'a' approaches zero.Hence, Equation ( 7) becomes (8) where p is the normalized horizontal coordinate of the model.
The trial solution of Equation ( 8) is u = −z(u;v)≈constant when ∂ 2 u/∂x 2 = 0 and ∂ 2 u/∂p 2 = 0.This solution means that all sliders are moving uniformly.Consider to add a small perturbation to u(p,x), that is, u(p,x) = −z(u;v) + u o exp(iqp + Xx), where u o is a very small quantity, q is the wave-number, and X = ~/~o is the dimensionless angular frequency.This gives ∂u/∂x = Xu o exp(iqp+Xx), ∂ 2 u/∂x 2 = X 2 u o exp(iqp+Xx), and ∂ 2 u/∂p 2 = −q 2 u o exp(iqp+Xx).Substituting those quantities into Equation ( 8) leads to This shows that any small perturbation added to a slider, no matter how long its wavelength is, does not diverge when u = −z(u;v)≈constant.

Slip-weakening friction
As mentioned above, the simplified slip-weaken- where D is the dimensionless characteristic slip distance and also the dimensionless decreasing rate of friction with slip.This makes Equation (8) become (9) In order to solve Equation (9), we take the Laplace Transformation (LT, denoted by L), which is described in Appendix A. According to Appendix A, the LT of Equation ( 9) is (10) This gives (11) To solve Equation (11), let U = U c +U p where U c and U p are, respectively, the complementary and particular solutions.Let } = (s 2 +1−D -1 ) 1/2 .According to the method shown in Appendix B, the general solution of Equation ( 11) is (12) where U c = C 1 e -}p/h +C 2 e }p/h and U p =−1/s} 2 .
Equation ( 11) consists of two kinds of waves: the first one is the travelling wave along the model represented by the first and second terms in its right-handed-side and the second one is the displacement at a site given by the third term.The first and second terms represent the waves travelling along the direction of increasing p and that of decreasing p, respectively.The second term with p< 0 can indeed be re-written as e -}|p|/h .From Equation (A2) in Appendix A, L -1 [e -}|p|/h ] with |p|/h>0 is (13 where C=C 1 or C 2 , c = (1−D -1 ) 1/2 and H(x−|p|/h) is the unit step function (H(x)=0 as x< 0 and H(x)=1 as x ≥ 0) representing a travelling plane wave.
The two parameters x =~ot and p= x/D o are, respectively, the normalized time and normalized (dimensionless) rupture distance, Hence, h is the normalized (dimensionless) rupture velocity and equal to v R /D o ~o, where v R is the rupture velocity.The arrival time of rupture wave at p or x is x a =|p|/h or t a =x/v R .When the rupture propagates from 0 to p L , which is the normalized (dimensionless) total rupture length and equal to L/D o (L = the total rupture length), the normalized total duration time is x D = p L /h, and thus the total duration time is Based on the above-mentioned parameters, Equation (13) can be transferred to (14) Hence, the FT of the first part of Equation ( 14), i.e., x c1 (y,t The FT of the second part of Equation ( 14), i.e., )]e -i~t dt, with an integral range from −∞ to +∞.At each site, the rupture occurs just between t a and t R = t a +T R where T R is the rise time of rupture.When t >T R .the rupture stops.We are only interested on the spectrum at a site and thus different values of t a at different sites just lead to phase changes and cannot influence the spectral amplitudes when the wave attenuation is ignored.Hence, the value of t a inside the integral can be set to be zero to simplify the problem.This gives X c2 (y,~)=−CD o c~ot a ∫[J 1 (c~ot)/t]e -i~t dt, with an integral range from 0 to T R .Since F [J 1 (t)/t] = [2(1 − ~2)/r] 1/2 rect(~/2c~o), where rect(~) is the rectangular function: rect(~) = 1 for |~|<1/2, 1/2 for |~|=1/2, and 0 for |~|>1/2, we have F [J 1 (c~ot)/t]= {2[1−(~/ c~o) 2 ]/r} 1/2 rect(~/2c~o).Clearly, the Fourier transform is equal to 1 only when |~|<c~o.This means that the Fourier amplitudes of x c2 (y,t) exist only when |~|<c~o.Hence, the complementary solution cannot lead to ~-n source scaling at high frequencies.
The third term of Equation ( 12), i.e., U p (p,s)= −1/2s(s 2 +1−D -1 ), represents source behavior of slip at a site.According to Appendix A, u p (p,x)=L -1 [U p (p,s)]= [1−cos(cx)]/2c.Obviously, the solution does not exist when 1−D -1 = 0 or D =1.When 1−D -1 ≠ 0 or D ≠ 1, we have (15) There are two cases in Equation ( 15): one with 1−D -1 <0 or D <1 and the other with 1−D -1 >0 or D >1.For the first case, sin The solution displays a negative hyperbolic sinetype function.It cannot represent a commonly-defined wave.For the second case, the solution Equation ( 14) shows vibrations, specified with a positive sine-type function, at a site.Replacing u p (p,x) by x p (y,t), p by x/D o , and x by ~ox into Equation (15) leads to (16) Obviously, x p (y,t) =0 when t = 0.The FT of x p (y,t), i.e., X p (y,~)=F [x p (y,t)], is (17) Obviously, the values of X p (y,~) exist only at ~= −c~o, 0, and c~o.thus unable to represent a commonly acceptable earthquake source displacement spectrum.This means that the simplified linearly slipweakening friction law cannot work well for earthquake ruptures based on the 1D BK model.Since this friction law has been widely used by others for other models as mentioned above, I assume that a more complicated slip-dependent friction law like that used by Cao and Aki [1984/85] or those proposed by Rice [2006] should be taken into account for the 1D BK model in the future.Bizzarri [2012] considered a single spring-slider system and assumed a linear slip-weakening friction law to govern the motion of the slider.He provided a closed-form, analytical solution of the physical system.His solution of slip includes two exponential time functions, like e -Ct for t>0.The function can be represented I by e -Ct H(t).Like the previous complimentary solution, the Fourier spectrum of his slip function is 1/(i~+C), thus showing ~-1 scaling in the entire range of ~.

Velocity-weakening friction
As mentioned above, the simplified velocity-weakening friction law is z(v)=1−v/(V c /D o ~o).Define the dimensionless characteristic velocity to be y =V c /D o ~o, we have z(v)=1−v/y, where y is also the dimensionless decreasing rate of friction with velocity.This makes Equation ( 8) become ( 18) According to Appendix A, the LT of Equation ( 18) is (19) This gives (20) Let g=(s 2 +1−sy -1 ) 1/2 .From Appendix B, the general solution of Equation ( 20) is ( 21) Equation ( 21) shows the superposition of two kinds of waves, i.e., the travelling waves along the model represented by the first and second terms in its right-handed-side and the displacement at a site shown by the third term.The first and second terms represent the waves travelling along the direction of increasing p and that of decreasing p, respectively.The second term with p < 0 can indeed be re-written as e -g|p|/h .Let Obviously, this expression represents a travelling plane wave.In the right-handed side, the quantities in the front of H(x−|p|/h) show the wave amplitude.Unlike Equation ( 13) for slip-weakening friction, Equation ( 22) has an additional term e x/2y , which increases exponentially with time as displayed in Figure 2. When the rupture propagates from 0 to p L , the normalized total rupture duration time is x D = p L /h, and thus the total rupture duration time is x=~ot, t a =y/v R , and s = ~o/2y into Equation ( 22) leads to (23) The first part of Equation ( 23) is (24) Its FT is X c1 (y,~)=F [x c1 (y,t)]=CD o ∫e -(i~-s)t dt.At each site, the integration is made from t a to t R =t a +T R , where T R is the rise time of rupture.Hence, we have The FT of the second part of Equation ( 23), i.e., )]e -i~t dt, with an integral range from −∞ to +∞.At each site, the rupture occurs just between t a and t R =t a +T R , where T R is the rise time of rupture.We are only interested on the spectrum at a site and thus different values of t a at different sites just lead to phase changes and cannot influence the spectral amplitudes when the wave attenuation is ignored.Hence, the value of t a inside the integral can be set to be zero to simplify the problem.This gives with an integral range from 0 to T R .Since F [ J 1 (t)/t] = [2(1−~2)/r] 1/2 rect(~/2), where rect(~) is the rectangular function as mentioned above, we have Obviously, the Fourier transform is equal to 1 when |~|<v~o.The Fourier transform of e st is F [e st ] = ∫e st e -i~t dt, with an integral range from 0 to T R , and equal to [exp(s −i~)T R −1]/(s −i~).Mathematically, F {[ J 1 (v~ot)/t]e st }=F [ J 1 (v~ot)/t] * F [e st ], where the symbol " * " denotes the convolution of F [ J 1 (v~ot)/t] and F [e st ], i.e., (26) with an integral range from −v~o to +v~o.Equation ( 26) is quite complicated and cannot be easily integrated.Since only spectral scaling at high frequencies is considered here, the integration is conducted only for high frequencies.Since |z|< v~o<~o and s =~oy/2<~o, s−i(~−z) can be approximated by i~when ~> >~o.Under this condition, Equation ( 26) can be re-written as ( 27) with an integral range from −v~o to +v~o.The two integrals in Equation ( 27) are finite and represented by I 1 and I 2 , respectively, for the former and latter.Hence, F {[ J 1 (v~ot)/t]e st }=(1/i~v~o)(2/r) 1/2 (I 1 e -i~− I 1 ).This gives e st }| is dependent upon ~, but it decreases with increasing ~in a power-law form of ~-1 .This means that the complementary solution shows ~-1 source scaling at high frequencies.Of course, the magnitude of |F {[ J 1 (v~ot)/t]e st }| is also in terms of v and ~o.
(ii) For 1−4y 2 <0 or y >0.5:This gives 4y 2 −1>0.Let q=(4y 2 −1) 1/2 , Consider a right triangle with three sides: the longest side with a length of L=(1 2 +q 2 ) 1/2 =2y, and the other two with lengths of q and 1, respectively.The angle between the longest side and the side with a length of 1 is set to be i.Hence, we have cos(i)=1/2y, sin(i)=q/2y, and tan(i)= q.The tangent term gives i = tan -1 (q).Define sin(a) = (e iqx/2y − e -iqx/2y )/2i and cos(a)=(e iqx/2y + e -iqx/2y )/2.Replacing u p (p,x) by J " e i i d , = -+ - Due to the existence of e st , the FT of the second part of Equation ( 30) cannot be obtained directly from the FT table.Hence, we must conduct the integration to get its FT.The integration is made in the range from 0 to T R .The second term of Equation ( 30) includes a sine function with a period of T p = 4ry/q~o.As displayed in Figure 4, the solid plus dotted curve represents the velocity waveform, while the solid curve denotes the displacement waveform.When t = T R , the displacement at a site reaches the maximum value.Since the wave is considered to propagate from −y to +y, the motion stops at the time instant when the velocity is equal to zero.In Figure 4, T p and T R are, respectively, the period of velocity waveform and the rise time of displacement waveform.Clearly, T R is roughly equal to T p /2, thus giving ~oT R /2y = r/q.The FT of e st sin(qst −i)/sini is made from 0 to T R , and the result is (31) The absolute value of X p2 (y,~) is (32) At ~= 0, |X p2 (y,~)|=(e 2rq +1)/~oy = r/~oy(4y 2 − 1) 1/2 , which is a function of both ~o and y.When (~o/y) 2 <<2, {~4+[(~o/y) 2 − 2]~2+~o 4 } 1/2 ≈~2−~o 2 .This leads to the existence of singularity when ~= ~o.Numerical tests show that singularity appears when ~o is lower than a certain angular frequency, ~l, which depends on y.Hence, in order to obtain a small value of ~l= ~o/y, the choice of ~o must be dependent on y.The larger y is, the higher ~o can be.For example, ~l= 0.4r when y =0.55.Because of ~o=(l l /t A ) 1/2 =(K l /m) 1/2 , weak coupling between the moving plate and the fault and/or a low density of fault-zone materials will lead to low ~o, thus resulting in singularity.
Essentially, Figure 5a is similar to Figure 5b.The spectral amplitude decreases with increasing y as well as increasing ~o.Higher y as well as larger ~o leads to smaller spectral amplitude.The difference in spectral amplitudes between two consequent values of y decreases with increasing y as well as increasing ~o.This means that at high y an increase in y cannot produce stronger spectral amplitudes.This phenomenon also exists for ~o.In Figures 5a and 5b There are some differences between Figure 5a and Figure 5b.First, the lines with a slope value of −1 at higher ~o are parallel to one another in Figure 5a and merge together at much higher ~o in Figure 5b.Secondly, ~c is almost independent on y (see Figure 5a) and increases with ~o (see Figure 5b).
A comparison between Equation ( 16) and Equation (29) reveals that a major difference between the two equations is the existence of an exponential term in Equation ( 29) and not in Equation ( 16).Such an additional effect is frequency-dependent.This makes the , e e i i .difference in the results between the two friction laws.Nonlinearity of slip-weakening friction, as shown in Equations ( 3)-( 5), could lead to different results.It is significant to examine the effect on source displacement spectrum caused by nonlinear slip-weakening friction in the future.
Figure 5 clearly show that linearly velocity-weakening friction results in |X p2 (y,~)|~~0 at low ~and |X p2 (y,~)|~~-1 at high ~.This is different from ~-square source scaling as proposed by Aki [1967].This might be due to a use of 1-D model in this study and a use of 2-D model by Aki [1967].The conventional concept is that ~-square scaling is caused by the smoothing effect due to time and fault length and ~-cubic scaling due to time, fault length, and fault width [cf.Aki and Richards 1980].However, the present model consists only of time and fault length.It seems that time does not causes the smoothing effect.I assume that ~-square scaling should be caused by the smoothing effect due to both fault length and fault width.This must be examined by using the 2-D model.
The present results are somewhat different from those obtained by Shaw [1993], even though using the same 1-D BK model in the two studies.Unlike his results, in this study a single ~-1 source scaling law at high ~is operative for both small and large earthquakes.A possible reason to cause the difference is that his velocity-weakening law is nonlinear and more complicated than the present one.Nonlinearity can produce unexpected results.It is significant to analytically and numerically investigate the effect of different friction laws on source scaling near future.
The previous studies as mentioned in the section of "Introduction" suggest that fractal and/or heterogeneous spatial distributions of fault strengths, displacements etc. on the fault plane are major factors in affecting scaling of earthquake source displacement spectrum.The present study obviously suggests that linearly velocity-weakening friction can also play a significant role on such scaling.The corner angular frequency (~c) is commonly considered to be a ratio of the rupture velocity to the fault length [cf.Aki and Richards 1980].Figure 5 exhibit that ~c increases with ~o, yet is not a function of y.This indicates that ~o, which depends on the mass of a slider and coupling between a slider and the moving plate, controls ~c, while the characteristic velocity or the decreasing rate of friction with velocity is not a factor in influencing ~c.
However, two problems should be resolved in the near future.First, the present study cannot suggest whether the ~-cubic source scaling holds or not, because only the 1-D model is taken into account.Secondly, the present result is valid only for the longitudinal-wave-type rupture.Included also in earthquake are the shear-wavetype waves.In order to explore the two problems, a similar study on the basis of the 2-D dynamical spring-slider mode extended from the 1-D BK model by Wang [2000Wang [ , 2012] ] should be performed further.

Conclusions
Scaling of earthquake source displacement spectra is analytically studied based on the continuous form of the one-dimensional dynamical spring-slider model proposed by Burridge and Knopoff [1967] in the pres- ence of two types of friction, i.e., linearly slip-weakening friction and linearly velocity-weakening friction.The general solution has the complementary and particular parts.The complementary solution exhibits ~-1 source displacement spectra for the two types of friction.For the particular solution, slip-weakening friction cannot produce ~-n source displacement spectra, while velocity-weakening friction results in acceptable spectra.For velocity-weakening friction, the log-log plot of spectral amplitude versus angular frequency is mainly controlled by the natural angular frequency, ~o, of the system and the (dimensionless) decreasing rate, y, of the friction law.Essentially, the spectral amplitude exhibits almost ~0 scaling at low ~and decreases with increasing ~following a line with a slope value of −1, which is the scaling exponent.However, the plots in the range of medium ~are clearly different between ~o<~l and ~o> ~l, where ~l depends on y.For example, ~o<0.4r and ~o> 0.4r when y= 0.55.For ~o> 0.4r, the spectral amplitude exhibits almost ~0 scaling when ĩs lower than the corner angular frequency (denoted by ~c), and the spectral amplitude monotonically decreases with ~and follows a line with a slope value of −1 when ~>~c.The corner angular frequency increases with ~o and independent on y.On the other hand, for ~o<0.4r, the spectral amplitude first exhibits almost ~0 scaling at low ~, then increases with ~up to a peak value, and finally decreases with increasing ~, following a line with a slope value of −1, which is the scaling exponent.The angular frequency associated with the peak value is almost ~o.Consequently, the source displacement spectrum shows ~-1 scaling at high ~when linearly velocity-weakening friction is taken into account under the condition of y >0.5.

Figure 2 .
Figure 2. The function of exp(~ot/2y) versus t.The solid curve from t a =y/v R to t R = t a +T R , where v R and T R are the rupture velocity and rise time, respectively, displays the integration range.
x p (y,t)/D o , p by y/D o , x by ~ot, and s = ~o/2y into Equation (29) leads to (30) Obviously, x p (y,t)=0 at t=0.This solution shows a sine-function-type standing wave.Unlike Equation (21), Equation (30) has an extra term e st , which shows a temporal increase in the amplitude with a constant of s -1 =2y/~o.Since the motion starts at t=0, x p (y,t) is null when t <0 and thus Equation (30) must have the form of x p (y,t)H(t).The FT of −D o H(t) is X p1 (y,~)=−D o [rd(~)− i/~].The first term of X p1 (y,~) is a delta function of slip at ~= 0, and it does not influence the source displacement spectrum when ~> 0. Since the second term of X p1 (y,~) exhibits a decay of motion with ~in a power-law function of ~-1 for the entire range of ~.
14r (from top to bottom) when y = 0.55.To produce Figure 5, D o is taken to be 1 unit.Since a uniform fault is considered here, the spectral amplitude is only a function of ~and related source parameters and independent on the position.Hence, |X (~)|=|X p2 (y,~)| is displayed in the figure.The values of |X (~)| are normalized by respective maximum values.The dashed line in each diagram has a slope value of −1. Figure , log(|X (~)|) is almost constant (exhibiting ~0 scaling) at low ~when ~is lower than a turning one, i.e., the corner angular frequency, ~c, which is located at the tip of the convex curve and not displayed in Figure 5, When ~> ~c, log(|X (~)|) monotonically decreases with ~following the respective lines with a slope value of −1, which is the scaling exponent.

Figure 4 .
Figure 4.The solid plus dotted curve represents the velocity waveform and the solid curve denotes the displacement waveform.T p and T R are, respectively, the predominant period of the velocity waveform and the rise time of displacement waveform.