A dynamical study of frictional effect on scaling of earthquake source displacement spectra ANNALS OF GEOPHYSICS, 59, 2, 2016, S0210; doi:10.4401/ag-6974 S0210 A dynamical study of frictional effect on scaling of earthquake source displacement spectra Jeen-Hwa Wang Institute of Earth Sciences, Academia Sinica, Nangang, Taipei, Taiwan ABSTRACT 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 fric- tion or linearly velocity-weakening friction. The main parameters of the model are the natural angular frequency, ~o, and the (dimensionless) de- creasing 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-weak- ening friction. The analytic solution includes the complementary and par- ticular 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 spec- tral 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- ponent. 1. Introduction The body-wave seismic spectrum, P(~), where ~= 2rf is the angular frequency and f is the frequency, is controlled by the seismic moment Mo 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 ap- propriate 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 displace- ment spectra from the seismograms recorded at nine near-fault stations generated by the 1999 Chi-Chi, Tai- wan, earthquake. Results show that the values of cor- ner frequency fc 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 scal- ing) 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 near- fault 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. [1987, 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 dimen- sion. In his model, an earthquake is composed of a hi- erarchical set of smaller events. The static stress drop is parameterized to rn, and strength is assumed to be pro- portional to static stress drop. He found that a distribu- tion 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 rup- ture area of the earthquake. Based on an ideal system under external random forces, Koyama and Hara [1993] Article history Received January 26, 2016; accepted April 8, 2016. Subject classification: Source spectrum, Spring-slider model, Weakening friction, Natural frequency. 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 expo- nent c is a function of the fractal dimension of the scal- ing 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 ve- locity 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 (abbre- viated by the 1-D BK model hereafter) proposed by Burridge and Knopoff [1967]. Meanwhile, friction con- trols 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 fric- tion, 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 Mo> Moc, where Mo is the seismic moment and Moc is the characteristic seismic moment, there are different power-law relations in three angular-frequency regions: P(~)~~0 for ~ < 2r/Mo; P(~)~~-1 for 2r/Mo< ~ <2r/p; and P(~)~~ -c for ~> 2r/p, where p= 2ln(4l2/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. Obvi- ously, there are two turning points in the source spec- trum 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 re- sult is somewhat between ~-2 and ~-3 models, because c is about 2.5. This is somewhat different from that pro- posed by Aki [1967]. For large events, Shaw’s theoreti- cal source spectra are more complicated than Aki’s. For small events with Mo< Moc, there are two power-law re- lations 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 theoret- ical result for small events is similar to Aki’s ~-square model. Shaw [1993] only used a velocity-weakening fric- tion law for numerical simulations. It seems necessary to explore the effect on the source displacement spec- trum due to slip-weakening friction. In this work, I will focus on the analytical study of frictional effect on scal- ing of earthquake source displacement spectra using the continuous form of the 1-D BK model in the pres- ence of two types of friction, i.e., slip-weakening fric- tion and velocity-weakening friction. 2. One-dimensional spring-slider model The 1-D BK model (see Figure 1) consists of N slid- ers of equal mass, m, and springs with one slider being linked by a coil spring of strength, Kc, with the other. WANG 2 Figure 1. An N-degree-of-freedom one-dimensional dynamical spring-slider system. 3 Each slider is also pulled by a leaf spring of strength, Kl, on a moving plate with a constant velocity, Vp. At time t = 0, all the sliders rest in the individual equilib- rium states. The i-th slider (i = 1, …, N) is located at po- sition Xi, measured from its initial equilibrium position, along the horizontal axis. Each slider is subjected to a slip- and/or velocity-dependent frictional force. The fric- tional force is denoted by Fi (Xi; Vi), where Vi=∂Xi/∂t, with a static frictional force, Fsi, at rest. Elastic strain of each slider gradually accumulates due to the moving plate. Once the elastic force at a slider is greater than Fsi, the static frictional force drops to the dynamic fric- tional force, Fdi, and the slider will move subject to Fdi. The equation of motion is (1) Wang [1995] defined the ratio l= Kc/Kl to be the stiffness ratio of the system. It is an important param- eter 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 Kl 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 Vp 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 lon- gitudinal- wave-type rupture [Wang 1996]. However, the result obtained from this model is still significant for un- derstanding the P-wave-type source scaling. In crack me- chanics, 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 consid- ered and therefore the problem can be treated easily. The present advantage is that the closed-form analytical solu- tions can be obtained. 3. Friction Velocity-dependent friction and slip-dependent friction The frictional force between two contact planes is classically considered to drop from static one to dy- namic 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. Al- though 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, weaken- ing-hardening friction law. Carlson and Langer [1989] and Carlson et al. [1991] considered a velocity-weaken- ing friction law: F(v) = Fo/(1 + v/vc), where Fo is the static frictional force and vc is a characteristic speed to specify the variation in F with velocity. F(v) is Fo at v = 0 and decreases monotonically from Fo to zero as |v| be- comes 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 mod- els are very unstable at low values of the fault slip ve- locity 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 scal- ing of source displacement spectra of earthquakes. An- alytic results will help us to understand such an effect. For the first-order approximation, Wang [1995, 1996] considered a piece-wise, linearly velocity-depen- dent weakening-hardening friction which is simplified from the friction law proposed by Burridge and Knopoff [1967]. The decreasing (weakening) and in- creasing (hardening) rates of dynamic friction strength with sliding velocity are two main parameters of this friction law. The two rates are denoted, respectively, by rw and rh. The function F(v) is defined only for v ≥ 0 and F(v) is a negative infinity when v < 0. This means that ; . m X X X X X X t K Kl Vp t Fi Vi 2i i i i i i Vi c Kl Vp Fi 2 2 1 12 2 = - + - - - + -R Q Q Q W V V V A DYNAMICAL STUDY OF EARTHQUAKE SPECTRA no backward motion is allowed. When v = vc, F(v) reaches the minimum value, i.e., gFo (0 0, X is purely imaginary. Hence, the solution is u(p,x)=−z(u;v) +uoexp[i(qp+Xx)]. 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- ing friction law is z(u)=1−u/(Xc/Do). Define D = Xc/Do, we have z(u)=1−u/D, where D is the dimensionless characteristic slip distance and also the dimensionless decreasing rate of friction with slip. This makes Equa- tion (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 = Uc+ Up where Uc and Up are, respectively, the complementary and particular solutions. Let } = (s2 + 1 − D-1)1/2. According to the method shown in Appendix B, the general solution of Equation (11) is (12) where Uc= C1e -}p/h + C2e }p/h and Up=−1/s} 2. Equation (11) consists of two kinds of waves: the first . m yd z yd z x x Kl yd z x yd z x t K x Fi 2 i i i i c i Kl Fi 2 2 1 12 2yd d yd d yd d yd d = - + - - + -Q R R R QV W W W V .x x x k xt k x fi2i i i i fic i lA 2 2 1 12 2t -= - + -+ -Q V ; .u u u u uu v2i i i i ii i 2 2 1 12 2x l z= - + - -+ -Q QV V ; .u a u u u a u u v2i i i i 1 i i i 2 2 2 1 22 2x l z= - + - --+Q QV V" % ; ,Vu u u u vh2 2 2 222 2 2 2x p z= - -R QW ,V .u h u u u12 2 2 2 22 2 2 2x p D= - - -R QW V , , , ,/ . s U s U s U s U s h s1 2 2 2 22 2p p p p p D = - - + Q Q Q Q V V V V / .h sU s U 112 2 2 2 12 2p D+ + - =-Q V , C eC e / ,U s s1hh 21 2 p }= + -}p }p-Q V A DYNAMICAL STUDY OF EARTHQUAKE SPECTRA 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 trav- elling along the direction of increasing p and that of de- creasing 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 = C1 or C2, 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/Do are, re- spectively, the normalized time and normalized (dimen- sionless) rupture distance, Hence, h is the normalized (dimensionless) rupture velocity and equal to vR/Do~o, where vR is the rupture velocity. The arrival time of rupture wave at p or x is xa=|p|/h or ta= x/vR. When the rupture propagates from 0 to pL, which is the nor- malized (dimensionless) total rupture length and equal to L/Do (L = the total rupture length), the normalized total duration time is xD= pL/h, and thus the total du- ration time is TD= xD/~o= pL/~oh, thus giving TD= pL/~oh = (L/Do)/~o(vR/Do~o) = L/vR. Based on the above-mentioned parameters, Equa- tion (13) can be transferred to (14) Hence, the FT of the first part of Equation (14), i.e., xc1(y,t) = CDoH (~ot−|y|/Doh), is Xc1(y,~)= F [xc1(y,t)]=F [CDoH(~ot−|y|/Doh)]=CDo[rd(~)− (i/~)exp(|y|/Doh)]. When ~ > 0, the spectrum shows ~-1 scaling. The FT of the second part of Equation (14), i.e., xc2(y,t) = −CDoc~otaJ1[c~o(t 2 − ta 2)1/2]/(t2 − ta 2)1/2} H[~o (t − ta)], is Xc2(y,~)=F [xc2(y,t)]=−CDoc~ota∫{ J1[c~o (t2 − ta 2)1/2]/(t2 − ta 2)1/2} H [~o(t − ta)]e -i~tdt, with an inte- gral range from −∞ to +∞. At each site, the rupture occurs just between ta and tR= ta+TR where TR is the rise time of rupture. When t >TR. the rupture stops. We are only interested on the spectrum at a site and thus differ- ent values of ta at different sites just lead to phase changes and cannot influence the spectral amplitudes when the wave attenuation is ignored. Hence, the value of ta inside the integral can be set to be zero to simplify the problem. This gives Xc2(y,~)=−CDoc~ota∫[J1(c~ot)/t]e -i~tdt, with an integral range from 0 to TR. Since F [J1(t)/t] = [2(1 − ~2)/r]1/2rect(~/2c~o), where rect(~) is the rectangu- lar function: rect(~) = 1 for |~|< 1/2, 1/2 for |~|= 1/2, and 0 for |~|>1/2, we have F [J1(c~ot)/t]= {2[1−(~/ c~o) 2]/r}1/2rect(~/2c~o). Clearly, the Fourier trans- form is equal to 1 only when |~| 0 or D >1. For the first case, sin [(1 −D-1)1/2 x/2] = i·sinh [(D-1−1)1/2 x/2], and thus sin2 [(1 −D-1)1/2 x/2] = − sinh2 [(D-1−1)1/2 x/2]. The solution displays a negative hyperbolic sine- type 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 up(p,x) by xp(y,t), p by x/Do, and x by ~ox into Equation (15) leads to (16) Obviously, xp(y,t) = 0 when t = 0. The FT of xp(y,t), i.e., Xp(y,~) =F [xp(y,t)], is (17) Obviously, the values of Xp(y,~) exist only at ~ = −c~o, 0, and c~o. thus unable to represent a com- monly acceptable earthquake source displacement spectrum. This means that the simplified linearly slip- weakening friction law cannot work well for earth- quake 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 com- plicated 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 sys- tem 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 func- tions, like e-Ct for t > 0. The function can be represented , C J H ,$ x y t t t t t t t t1c o a o o a a a 2 2 1V 2 1 2 2 1V 2c~ c~ ~- - = - -Q Q Q QV 1V V 1V ! " ,$ %E I , C JV H ,V u h h h h 1c 1 2 2V 1V 2 2 2V 1$ 2 p x c p c x p x p x p = - - - - Q Q Q Q QV 2 V JV ,V 2V 1V ! " 1$ %E I , .cos sinu 1 22p 2 2 p x cx c cx c-= =Q Q QV V V! $ , .sinu D t 2p o o 2 2 p x c~ c=Q RV W , . Xp y D 2 2 2 Xp o o o 2 ~ r d ~ d ~ c~ d ~ c~ c = - + + - Q Q Q QV V V V! $ E I WANG 6 7 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-weak- ening friction law is z(v) = 1− v/(Vc/Do~o). Define the dimensionless characteristic velocity to be y = Vc/Do~o, we have z(v) = 1− v/y, where y is also the dimension- less 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= (s2 + 1 − sy-1)1/2. From Appendix B, the gen- eral 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 v = (1 −1/4 y2)1/2 . According to Appendix A, L-1[e-g|p|/h] with |p|/h > 0 is (22) 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. Un- like Equation (13) for slip-weakening friction, Equation (22) has an additional term ex/2y, which increases ex- ponentially with time as displayed in Figure 2. When the rupture propagates from 0 to pL, the normalized total rupture duration time is xD= pL/h, and thus the total rupture duration time is TD= xD/~o= pL/~oh = (L/Do)/~oh. Substituting uc(p,x) = xc(y,t)/Do, p= y/Do, x=~ot, ta=y/vR, and s = ~o/2y into Equation (22) leads to (23) The first part of Equation (23) is (24) Its FT is Xc1(y,~)=F [xc1(y,t)]=CDo∫e -(i~ - s)tdt. At each site, the integration is made from ta to tR= ta+TR, where TR is the rise time of rupture. Hence, we have (25) This give |Xc1(y,~)|= CDo|exp(stR) − exp(sta)|/ (~2+s2)1/2. Clearly, the spectral amplitude decays with ~-1 in the entire range of ~, and |Xc1(y,0)|= CDo|exp (stR) − exp(sta)|/s. The FT of the second part of Equation (23), i.e., xc2(y,t) = −CDov~otae st J1[v~o(t 2 − ta 2)1/2] H [~o(t − ta)]/(t 2 − ta 2)1/2}, is Xc2(y,~) =F [xc2(y,t)] = −CDov~ota ∫{ J1[v~o(t 2 − ta 2)1/2]/(t2 − ta 2)1/2} est [H [~o(t − ta)] e -i~tdt, with an integral range from −∞ to +∞. At each site, the rupture occurs just between ta and tR=ta+TR, where TR is the rise time of rupture. We are only interested on the spectrum at a site and thus different values of ta at different sites just lead to phase changes and cannot influence the spectral amplitudes when the wave at- tenuation is ignored. Hence, the value of ta inside the integral can be set to be zero to simplify the problem. This gives Xc2(y,~) = −CDov~ota∫[ J1(v~ot)/t] e st e-i~tdt, .u h u u v12 2 2 2 22 2 2 2x p y= - - -R QW V , , , ,/ s U s h U s U s U ss s1 2 2 2 22 2p p p p p y = - - +Q Q Q QV V V V /h U s U ss 112 2 2 2 12 2p y =- + - -Q V , C e C e / .U s s1h h1 2 2 p g= + -gp gp-Q V , C JV H . eu h h h h 1c 2 1 2 2V 1V 2 2 2V 1 2 p x v p v x p x p x p = - - - - x yQ Q Q Q Q QV 2 V JV V 2V 1V ! " $ %E I , C e JV H . x y t D t t t t t t t 1c o o a o a / a / o a t 1 2 2 1V 2 2 2 1V 2 v~ v~ ~ = - - - - sQ Q Q Q QV 1V JV V 1V ! " $ %E I , C e H .x y t D t tc o o a t 1 ~= - sQ QV V! $ , C . exp exp X y i i D i t t c o o o a R1 ~ ~ s ~ s ~ s- + - + = - + -Q Q Q Q V V V V! ! $ $E H A DYNAMICAL STUDY OF EARTHQUAKE SPECTRA Figure 2. The function of exp(~ot/2y) versus t. The solid curve from ta= y/vR to tR= ta+TR, where vR and TR are the rupture veloc- ity and rise time, respectively, displays the integration range. with an integral range from 0 to TR. Since F [ J1(t)/t] = [2(1− ~2)/r]1/2 rect(~/2), where rect(~) is the rectan- gular function as mentioned above, we have F [ J1(v~ot)/t] ={2[1−(~/v~o) 2]/r}1/2 rect(~/2v~o). Obviously, the Fourier transform is equal to 1 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 I1 and I2, respectively, for the former and latter. Hence, F{[ J1(v~ot)/t] e st } = (1/i~v~o)(2/r) 1/2(I1e -i~− I1). This gives |F {[ J1(v~ot)/t] e st } |=|~v~o| -1(2/r)1/2|[(I1 2 − 1)cos2(~) − 2I1I2cos(~)+I2 2]1/2| Although the magnitude of |F {[ J1(v~ot)/t] 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 mag- nitude of |F {[ J1(v~ot)/t] e st } | is also in terms of v and ~o. The third term of Equation (21), i.e., Up(p,s) = −1/2s(s2+1− sy-1), directly shows source behavior of slip at a site and can be re-written as −(1−4y2)-1/2{−1/(s+ a)+[1−(1−4y2)1/2]/2(s+b)−[1+(1−4y2)1/2]/2(s+c)}, where a = 0, b = − (1/2y)[1+(1−4y2)1/2], and c= − (1/2y) [1−(1−4y2)1/2]. Obviously, the solution does not exists when 1−4y2 = 0 or y= 0.5. When 1−4y2≠0 or y≠0.5, the ILT of Up(p,s) is In Equation (28), there are two cases: one with 1−4y2 >0 or y<0.5 and the other with 1−4y2 <0 or y> 0.5. The respective results are described below. (i) For 1−4y2 >0 or y<0.5: Let sinh(a) = (ecx/2y − e-cx/2y)/2 and cosh(a) = (ecx/2y +e-cx/2y)/2. Equation (28) is re-written as (29) Obviously, the solution cannot represent a com- monly known wave and thus is excluded in this study. (ii) For 1−4y2 <0 or y> 0.5: This gives 4y2 −1> 0. Let q = (4y2 −1)1/2, Consider a right triangle with three sides: the longest side with a length of L = (12 + q2)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) = (eiqx/2y − e-iqx/2y)/2i and cos(a)=(eiqx/2y + e-iqx/2y)/2. Replacing up(p,x) by xp(y,t)/Do, p by y/Do, x by ~ot, and s = ~o/2y into Equation (29) leads to (30) Obviously, xp(y,t) = 0 at t = 0. This solution shows a sine-function-type standing wave. Unlike Equation (21), Equation (30) has an extra term est, which shows a temporal increase in the amplitude with a constant of s-1 = 2y/~o. Since the motion starts at t = 0, xp(y,t) is null when t < 0 and thus Equation (30) must have the form of xp(y,t)H(t). The FT of −DoH(t) is Xp1(y,~)=−Do[rd(~)− i/~]. The first term of Xp1(y,~) is a delta function of slip at ~= 0, and it does not influence the source dis- placement spectrum when ~> 0. Since the second term of Xp1(y,~) exhibits a decay of motion with ~ in a power-law function of ~-1 for the entire range of ~. J" e i i d , F exp t t z z TR z 2 1 1 o o / / t TR 1 1 V 2 2W 1 2 # v~ r v~ s ~ s ~- = - - - - s Q Q R Q Q 1 V V V 2W V " " J" ! % % % $%F F I I# J! e i e i d d , F exp t t z TR z z z z 2 o o / i o / o / t TR 1 1 1V 2 2V 2 1 2 2V 2 1 2 v~ ~v~ r v~ v~ = - - - - s ~- -Q Q Q Q Q QV V 1 V V 2V 2V J! ! ! $ % $ E E H H # # , e e e e . eu 2 2p 2 2 2 2 2 p x c c c- =- + - +cx y cx y cx y cx y cx y - - Q Q QV V V" % , e .sinh coshu 2p 2 p x c a c a c=- + -x yQ Q QV V V! $F I , e .sin sinx y t q tD 1p o t s i i+ -=- sQ Q QV V V" % WANG 8 Figure 3. The plot of ~*/~o versus v.(28) 9 Meanwhile, its amplitude is only Do, it cannot result in a remarkable effect on the source displacement spectrum. Define ~*= q~o= (4y 2 −1)1/2~o be the predominant an- gular frequency of the wave. Figure 3 displays the plot of ~* versus ~o. Obviously, ~ * first increases with ~o at low ~o and then slightly increases with ~o at high ~o. Due to the existence of est, 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 TR. The second term of Equation (30) includes a sine function with a period of Tp= 4ry/q~o. As displayed in Figure 4, the solid plus dotted curve represents the ve- locity waveform, while the solid curve denotes the dis- placement waveform. When t = TR, 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, Tp and TR are, respectively, the period of velocity waveform and the rise time of displacement waveform. Clearly, TR is roughly equal to Tp/2, thus giving ~oTR/2y = r/q. The FT of e st sin (qst −i)/sini is made from 0 to TR, and the result is (31) The absolute value of Xp2(y,~) is (32) At ~= 0, |Xp2(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. Nu- merical tests show that singularity appears when ~o is lower than a certain angular frequency, ~l, which de- pends 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= (ll/tA) 1/2= (Kl/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. Figure 5 displays the log-log plot of |Xp2(~)| ver- sus ~: (a) for y = 0.55, 0.60, 0.65, 0.70, 0.75, 0.80, 0.85, 0.90, 0.95, and 1.00 (from top to bottom) when ~o= 2r; and (b) for ~o= 0.2r, 0.4r, r, 2r, 4r, 6r, 8r, 10r, 12r, and 14r (from top to bottom) when y = 0.55. To pro- duce Figure 5, Do 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 in- dependent on the position. Hence, |X (~)|=|Xp2(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 5b shows a bump around ~ = ~o when ~o=0.2r as men- tioned previously. 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 de- creases 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, log(|X (~)|) is almost constant (exhibiting ~0 scaling) at low ~ when ~ is lower than a turning one, i.e., the corner angular fre- quency, ~c, which is located at the tip of the convex curve and not displayed in Figure 5, When ~ > ~c, log(|X (~)|) monotonically decreases with ~ follow- ing the respective lines with a slope value of −1, which is the scaling exponent. 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. Sec- ondly, ~c is almost independent on y (see Figure 5a) and increases with ~o (see Figure 5b). A comparison between Equation (16) and Equa- tion (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 addi- tional effect is frequency- dependent. This makes the , e e i i . Xp cos sin y D TR TR 1 2 Xp o o o t TR t TR 2 2 2 2 ~ ~ ~ ~ y ~ ~ ~ s~ = + - - - + s s Q Q R Q Q V V W V V ! $ , e e . Xp cos y D TR 1 2 Xp o / o / o o / TR q q 2 2 $ 2 2W 2 1 2 4 2W 2 2 1 2 ~ ~ ~ y ~ ~ ~ y ~ ~ = + + + + - r r Q R Q RV 2W V 2W! " "1$ % % F I A DYNAMICAL STUDY OF EARTHQUAKE SPECTRA Figure 4. The solid plus dotted curve represents the velocity wave- form and the solid curve denotes the displacement waveform. Tp and TR are, respectively, the predominant period of the velocity waveform and the rise time of displacement waveform. 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-weak- ening friction results in |Xp2(y,~)|~~ 0 at low ~ and |Xp2(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 re- sults, 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 veloc- ity- weakening law is nonlinear and more complicated than the present one. Nonlinearity can produce unex- pected results. It is significant to analytically and nu- merically 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 hetero- geneous spatial distributions of fault strengths, dis- placements etc. on the fault plane are major factors in affecting scaling of earthquake source displacement spec- trum. 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 ve- locity 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-wave- type waves. In order to explore the two problems, a sim- ilar study on the basis of the 2-D dynamical spring-slider mode extended from the 1-D BK model by Wang [2000, 2012] should be performed further. 5. Conclusions Scaling of earthquake source displacement spec- tra 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- WANG 10 Figure 5. The log-log plot of |X(~)| versus ~: (a) for y= 0.55, 0.60, 0.65, 0.70, 0.75, 0.80, 0.85, 0.90, 0.95, and 1.00 (from top to bottom) when Do= 1 unit and ~o=2r; and (b) for ~o= 0.2r, 0.4r, r, 2r, 4r, 6r, 8r, 10r, 12r, and 14r (from top to bottom) (from top to bottom) when Do= 1 unit and y= 0.55. 11 ence of two types of friction, i.e., linearly slip-weaken- ing friction and linearly velocity-weakening friction. The general solution has the complementary and par- ticular parts. The complementary solution exhibits ~-1 source displacement spectra for the two types of fric- tion. For the particular solution, slip-weakening friction cannot produce ~-n source displacement spectra, while velocity-weakening friction results in acceptable spec- tra. 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 ex- hibits almost ~0 scaling at low ~ and decreases with in- creasing ~ 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 ~ is lower than the corner angular frequency (denoted by ~c), and the spectral amplitude monotonically de- creases with ~ and follows a line with a slope value of −1 when ~>~c. The corner angular frequency in- creases 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. Acknowledgements. The author thanks an anonymous re- viewer and Prof. Andrea Bizzarri (Associated Editor of the journal) for useful comments and suggestions to improve the article. This study was financially supported by Academia Sinica, the Ministry of Science and Technology (Grand Nos.: MOST 103-2116-M-001-010 and MOST 104-2116-M-001-007), and the Central Weather Bureau (Grand Nos.: MOTC-CWB-104-E-07 and MOTC-CWB-105-E-02), ROC. References Aki, K. (1967). Scaling law of seismic spectrum, J. Geo- phys. Res., 72, 1217-1231. Aki, K. (1972). Scaling law of earthquake source-time function, Geophys. J. R. Astro. Soc., 31, 3-25. Aki, K. (1979). Characterization of barriers on earth- quake fault, J. Geophys. Res., 84, 6140-6148. Aki, K., and P.G. Richards (1980). QUANTATIVE SEIS- MOLOGY: Theory and Methods, WH Freeman and Company, 932 p. Bak, P., C. Tang and K. Wiesenfeld (1987). Self-orga- nized criticality: An explanation of 1/f noise, Phys. Rev. Lett., 59, 381-384. Bak, P., C. Tang and K. Wiesenfeld (1988). Self-orga- nized criticality, Phys. Rev. A, 38, 364-374. Beeler, N.M., T.E. Tullis and D.L. Goldsby (2008). Con- stitutive relationships and physical basis of fault strength due to flash heating, J. Geophys. Res., 113, B01401; doi:10.1029/2007JB004988. Bernard, P., A. Herrero and C. Berge (1996). Modeling directivity of heterogeneous earthquake ruptures, B. Seismol. Soc. Am., 86, 1149-1160. Bizzarri, A., and M. Cocco (2006a). A thermal pressur- ization model for the spontaneous dynamic rupture propagation on a 3-D fault: 1. Methodological ap- proach, J. Geophys. Res., 111, B05303; doi:10.1029/ 2005JB003862. Bizzarri, A., and M. Cocco (2006b). A thermal pressur- ization model for the spontaneous dynamic rupture propagation on a 3-D fault: 2. Traction evolution and dynamic parameters, J. Geophys. Res., 111, B05304; doi:10.1029/2005JB003864. Bizzarri, A., and M. Cocco (2006c). Comment on “Earth- quake cycles and physical modeling of the process leading up to a large earthquake”, Earth Planet Space, 58, 1525-1528. Bizzarri, A. (2011a). On the deterministic description of earthquakes, Rev. Geophys., 49, RG3002; doi:10.1029/ 2011RG000356. Bizzarri, A. (2011b). Modeling repeated slip failures on faults governed by slip- weakening friction, B. Seis- mol. Soc. Am., 102, 812-821; doi:10.1785/0120110141. Bizzarri, A. (2012). What can physical source models tell us about the recurrence time of earthquakes?, Earth Sci. Rev., 115, 304-318; doi:10.1016/j.earscirev.2012. 10.004. Boatwright, J. (1978). Detailed spectral analysis of two small New York state earthquakes, B. Seismol. Soc. Am., 68, 1117-1131. Brune, J.N. (1970). Tectonic stress and the spectra of seis- mic shear waves from earthquakes, J. Geophys. Res., 75, 4997-5009. Burridge, R., and L. Knopoff (1967). Model and theo- retical seismicity. B. Seismol. Soc. Am., 57, 341-371. Cao, T., and K. Aki (1984/85). Seismicity simulation with a mass-spring model and a displacement hard- ening-softening friction law, Pure Appl. Geophys., 122, 10-24. Cao, T., and K. Aki (1986). Seismicity simulation with a rate- and state-dependent friction law, Pure Appl. Geophys., 124, 487-513. Carlson, J.M., and J.S. Langer (1989). Mechanical model of an earthquake fault, Phys. Rev. A, 40, 6470-6484. Carlson, J.M. (1991). Time intervals between charac- teristic earthquakes and correlation with smaller A DYNAMICAL STUDY OF EARTHQUAKE SPECTRA events: An analysis based on a mechanical model of a fault, J. Geophys. Res., 96, 4255-4267. Carlson, J.M., J.S. Langer, B.E. Shaw and C. and Tang (1991). Intrinsic properties of a Burridge- Knopoff model of an earthquake fault, Phys. Rev. A, 44, 884-897. Cochard, A., and R. Madariaga (1994). Dynamic fault- ing under rate-dependent friction, Pure Appl. Geo- phys., 142, 419-445; doi:10.1007/BF00876049. Das, S., and K. Aki (1977). Fault plane with barriers: a versatile earthquake model, J. Geophys. Res., 82, 5658-5670. Dieterich, J.H. (1972). Time-dependent friction in rocks, J. Geophys. Res., 77, 3690-3697. Dieterich, J.H. (1979). Modeling of rock friction, J. Geo- phys. Res., 84, 2161-2169. Dysart, P.S., J.A. Snoke and I.S. Sacks (1988). Source pa- rameters and scaling relations for small earthquakes in the Matsushiro region, Southwest Honshu, Japan, B. Seismol. Soc. Am., 78, 571-589. Fialko, Y.A. (2004). Temperature fields generated by the elastodynamic propagation of shear cracks in the Earth, J. Geophys. Res., 109, B01303; doi:10.1029/ 2003JB002496. Frankel, A. (1991). High-frequency spectral fall off of earthquakes, fractal dimension of complex rupture, b value, and the scaling of strength on faults, J. Geo- phys. Res., 96, 6291-6332. Haskell, N. (1966). Total energy and energy spectral density of elastic wave radiation from propagation faults, 2, A statistical source model, B. Seismol. Soc. Am., 56, 125-140. Herrero, A., and P. Bernard (1994). A kinematic self- similar rupture process for earthquakes, B. Seismol. Soc. Am., 84, 1216-1228. Hisada, Y. (2000). A theoretical omega-square model considering the spatial variation in slip and rupture velocity, B. Seismol. Soc. Am., 92, 387-400. Huang, M.W., and J.H. Wang (2002). Scaling of ampli- tude spectra of near-fault seismograms of the 1999 Chi-Chi, Taiwan, earthquake, Geophys. Res. Lett., 29 (10), 84-1 - 84-4. Johnson, R.E., and F.L. Kiokemeister (1968). Calculus with Analytic Geometry, Allyn and Bacon Inc., Boston, Mass, USA, 798 p. Kanamori, H. and D.L. Anderson (1975). Theoretical basis of some empirical relations in seismology, B. Seismol. Soc. Am., 65, 1073-1095. Kanamori, H., and G.S. Stewart (1978). Seismological aspects of Guatemala earthquake of February 4, 1976, J. Geophys. Res., 83, 3427-3434. Kim, S.G., Y.T. Chen, Z.L. Wu and G.I. Panza (1997). A mathematical theorem useful for the direct estima- tion of seismic source spectra, B. Seismol. Soc. Am., 87, 1281-1287 Knopoff, L., J.A. Landoni and M.S. Abinante (1992). Dy- namical model of an earthquake fault with local- ization, Phys. Rev. A, 46, 7445-7449. Koyama, J., and H. Hara (1993). Scaled Langevin equa- tion to describe the 1/fa spectrum, Phys. Rev. A, 46, 1844-1849. Madariaga, R., and A. Cochard (1994). Seismic source dynamics, heterogeneity and friction, Annali di Ge- ofisica, 37 (6), 1349-1375. Marone, C. (1998). Laboratory-derived friction laws and their application to seismic faulting, Annu. Rev. Earth Planet. Sci., 26, 643-696. Nur, A. (1978). Nonuniform friction as a physical basis for earthquake mechanics, Pure Appl. Geophys., 116, 964-989. Nussbaum, J., and A. Ruina (1987). A two degree-of-free- dom earthquake model with static/dynamic friction, Pure Appl. Geophys., 125, 629-656. Ohnaka, M. (2000). A physical scaling relation between the size of an earthquake and its nucleation zone size, Pure Appl. Geophys., 157, 2259-2282. Patanè, D., F. Ferrucci, E. Giampiccolo and L. Scara- muzzino (1997). Source scaling of microearthquakes at Mt. Etna volcano and in the Calabrian Arc (south- ern Italy), Geophys. Res. Lett., 24, 1879-1882. Rice, J.R. (1993). Spatio-temporal complexity of slip on a fault, J. Geophys. Res., 98, 9885-9907. Rice, J.R. (2006). Heating and weakening of faults dur- ing earthquake slip, J. Geophys. Res., 111, B05311; doi:10.1029/2005JB004006. Ruina, A.L. (1983). Slip instability and state variable fric- tion laws, J. Geophys. Res., 88, 10359-10370. Shaw, B.E. (1993). Moment spectra in a simple model of an earthquake fault, Geophys. Res. Lett., 20, 643-646. Shimamoto, T. (1986). Transition between frictional slip and ductile flow for halite shear zones at room tem- perature, Science, 231, 711-714. Wang, J.H. (1995). Effect of seismic coupling on the scaling of seismicity, Geophys. J. Int., 121, 475-488. Wang, J.H. (1996). Velocity-weakening friction law as a factor in controlling the frequency-magnitude rela- tion of earthquakes, B. Seismol. Soc. Am., 86, 701- 713. Wang, J.H. (1997). Effect of frictional healing on the scaling of seismicity, Geophys. Res. Lett., 24, 2527- 2530. Wang, J.H. (2000). Instability of a two-dimensional dy- namic spring-slider model of earthquake faults, Geophys. J. Int., 143, 389-394. Wang, J.H. (2002). A dynamical study of two one-state- variable, rate-dependent and state-dependent fric- WANG 12 13 tion laws, B. Seismol. Soc. Am., 92, 687-694. Wang, J.H. (2004). The seismic efficiency of the 1999 Chi-Chi, Taiwan, earthquake, Geophys. Res. Lett., 31, L10613; doi:10.1029/204GL019417. Wang, J.H. (2006a). Energy release and heat generation during the 1999 Chi-Chi, Taiwan, earthquake, J. Geo- phys. Res., 111, B11312; doi:10.1029/2005JB004018. Wang, J.H. (2006b). A review of the source parameters of the 1999 Chi-Chi, Taiwan, earthquake, Terr. Atmos. Ocean. Sci., 17, 179-202. Wang, J.H. (2007). A dynamic study of the frictional and viscous effects on earthquake rupture: a case study of the 1999 Chi-Chi earthquake, Taiwan, B. Seismol. Soc. Am., 97 (4), 1233-1244. Wang, J.H. (2008). One-dimensional dynamical model- ing of earthquakes: A review, Terr. Atmos. Ocean. Sci., 19, 183-203. Wang, J.H. (2009). Effect of thermal pressurization on ra- diation efficiency, B. Seismol. Soc. Am., 99, 2293-2304. Wang, J.H. (2012). Some intrinsic properties of the two- dimensional dynamical spring-slider model of earth- quake faults, B. Seismol. Soc. Am., 102, 822-835. Wang, J.H. (2013). Stability analysis of slip of a one- body spring-slider model in the presence of thermal pressurization, Annals of Geophysics, 56 (3), R03332; doi:10.4401/ag-5548. Wang, J.H. (2016). Studies of earthquake energies in Taiwan: A Review, Terr. Atmos. Ocean. Sci., 27 (1), 1- 19; doi:10.3319/TAO.2015.10.13.01(T). * Corresponding author: Jeen-Hwa Wang, Institute of Earth Sciences, Academia Sinica, Nangang, Taipei, Taiwan; email: jhwang@earth.sinica.edu.tw. © 2016 by the Istituto Nazionale di Geofisica e Vulcanologia. All rights reserved. A DYNAMICAL STUDY OF EARTHQUAKE SPECTRA Appendix A: The Laplace transformation The Laplace transformation (LT, denoted by L), of a function u, that is, L [u(p,x)] = U(p,s)=∫e-sxu(p,x)dx, in which the integration is made from –∞ to +∞. The Laplace transforms of several functions used in this study are: L [d2u(p,x)/dx2] = s2U(p,s), L [du(p,x)/dx] = sU(p,s), and L [1] = 1/s. In order to obtain the inverse Laplace transform (ILT, L-1) of U(s)=1/(s+a)(s+b)(s+c), which is the main form in this study, U(s) is first trans- ferred to −[(b−c)/(s+a)+(c−a)/(s+b)+(a−b)/(s+c)]/ (a−b)(b−c)(c−a). Due to L-1[1/(s+d)] = e-dx, (A1) Also used in this study is the following formula: (A2) where J1[y] and H [y] denote the Bessel function of the first kind with order 1 and the unit step function, respec- tively, when a>0, In addition, L-1{U(s −a)} = eaxu(x) and L-1[1/s] =1 are also used in this study. Appendix B: The solution of a second-order linear differential equation The method to solve a second-order linear differ- ential equation can be seen in Johnson and Kiokemeister [1968]. For a second-order inhomogeneous linear dif- ferential equation: (B1) where a, b, and c are three constants, the solution is U = Uc+ Up where Uc and UP are, respectively, the com- plementary and particular solutions. Uc is the solution of a homogeneous form of Equation (B1), i.e., (B2) Inserting the trial solution Uc= e mp into Equation (B2) leads to am2emp+ bmemp+ cemp= 0. This gives am2 + bm + c = 0. Thus, the solutions of m are m = −b± (b2−4ac)1/2/2a. Hence, the general form of Uc is C1e -mp+ C2e mp. Let y (p) = emp and z (p) = e-mp. To solve Up, we let (B3) The two function y1 and z1 should meet the fol- lowing conditions: y1’y + z1’z = 0 and y1’y ’+ z1’z’= f (x), where y1’=dy1/dp, y’=dy/dp, z1’=dz1/dp, and z’=dz/dp. Hence, y1 and z1 can be solved from the following in- tegral equations: y1= −∫[zf (x)(y1’y ’+ z1’z)/(yz’−y’z)]dp and z1=∫[yf (x)(y1’y ’+ z1’z)/(yz’−y’z)]dp. Inserting the related quantities into Equation (B3) leads to the par- ticular solution. b c e c a e a b e a b b c c a . L U s a b c t t t 1 =- - + - + - - - - - - - -Q Q Q Q Q Q Q V V V V V V V ! !$ $ Jb b H ,V L exp s 1 / / / 1 2 2 1V 2 1 2 2 1V 2 2 2 1V 2 a b a Jb x a x a x a - + = - - - - - Q Q Q Q 1V 1V 1V ,V! " $ % F F I I ad d dd c ,VbU x U x x U x f x2 2p + + =Q Q Q QV V V ,V ad d bd d c .U U U 02 2 +p p p p p+ =Q Q QV V V .Up y yV z zUp 1 1p p p p p= +Q Q Q Q QV yV V V V WANG 14 << /ASCII85EncodePages false /AllowTransparency false /AutoPositionEPSFiles false /AutoRotatePages /None /Binding /Left /CalGrayProfile (Dot Gain 20%) /CalRGBProfile (sRGB IEC61966-2.1) /CalCMYKProfile (U.S. Web Coated \050SWOP\051 v2) /sRGBProfile (sRGB IEC61966-2.1) /CannotEmbedFontPolicy /Warning /CompatibilityLevel 1.3 /CompressObjects /Tags /CompressPages true /ConvertImagesToIndexed true /PassThroughJPEGImages true /CreateJobTicket false /DefaultRenderingIntent /Default /DetectBlends true /DetectCurves 0.1000 /ColorConversionStrategy /LeaveColorUnchanged /DoThumbnails false /EmbedAllFonts true /EmbedOpenType false /ParseICCProfilesInComments true /EmbedJobOptions true /DSCReportingLevel 0 /EmitDSCWarnings false /EndPage -1 /ImageMemory 1048576 /LockDistillerParams true /MaxSubsetPct 100 /Optimize false /OPM 1 /ParseDSCComments true /ParseDSCCommentsForDocInfo true /PreserveCopyPage true /PreserveDICMYKValues true /PreserveEPSInfo true /PreserveFlatness true /PreserveHalftoneInfo false /PreserveOPIComments false /PreserveOverprintSettings true /StartPage 1 /SubsetFonts true /TransferFunctionInfo /Apply /UCRandBGInfo /Preserve /UsePrologue false /ColorSettingsFile (None) /AlwaysEmbed [ true /AndaleMono /Apple-Chancery /Arial-Black /Arial-BoldItalicMT /Arial-BoldMT /Arial-ItalicMT /ArialMT /CapitalsRegular /Charcoal /Chicago /ComicSansMS /ComicSansMS-Bold /Courier /Courier-Bold /CourierNewPS-BoldItalicMT /CourierNewPS-BoldMT /CourierNewPS-ItalicMT /CourierNewPSMT /GadgetRegular /Geneva /Georgia /Georgia-Bold /Georgia-BoldItalic /Georgia-Italic /Helvetica /Helvetica-Bold /HelveticaInserat-Roman /HoeflerText-Black /HoeflerText-BlackItalic /HoeflerText-Italic /HoeflerText-Ornaments /HoeflerText-Regular /Impact /Monaco /NewYork /Palatino-Bold /Palatino-BoldItalic /Palatino-Italic /Palatino-Roman /SandRegular /Skia-Regular /Symbol /TechnoRegular /TextileRegular /Times-Bold /Times-BoldItalic /Times-Italic /Times-Roman /TimesNewRomanPS-BoldItalicMT /TimesNewRomanPS-BoldMT /TimesNewRomanPS-ItalicMT /TimesNewRomanPSMT /Trebuchet-BoldItalic /TrebuchetMS /TrebuchetMS-Bold /TrebuchetMS-Italic /Verdana /Verdana-Bold /Verdana-BoldItalic /Verdana-Italic /Webdings ] /NeverEmbed [ true ] /AntiAliasColorImages false /CropColorImages true /ColorImageMinResolution 150 /ColorImageMinResolutionPolicy /OK /DownsampleColorImages true /ColorImageDownsampleType /Bicubic /ColorImageResolution 300 /ColorImageDepth -1 /ColorImageMinDownsampleDepth 1 /ColorImageDownsampleThreshold 1.10000 /EncodeColorImages true /ColorImageFilter /DCTEncode /AutoFilterColorImages true /ColorImageAutoFilterStrategy /JPEG /ColorACSImageDict << /QFactor 0.15 /HSamples [1 1 1 1] /VSamples [1 1 1 1] >> /ColorImageDict << /QFactor 0.15 /HSamples [1 1 1 1] /VSamples [1 1 1 1] >> /JPEG2000ColorACSImageDict << /TileWidth 256 /TileHeight 256 /Quality 30 >> /JPEG2000ColorImageDict << /TileWidth 256 /TileHeight 256 /Quality 30 >> /AntiAliasGrayImages false /CropGrayImages true /GrayImageMinResolution 150 /GrayImageMinResolutionPolicy /OK /DownsampleGrayImages true /GrayImageDownsampleType /Bicubic /GrayImageResolution 300 /GrayImageDepth -1 /GrayImageMinDownsampleDepth 2 /GrayImageDownsampleThreshold 1.10000 /EncodeGrayImages true /GrayImageFilter /DCTEncode /AutoFilterGrayImages true /GrayImageAutoFilterStrategy /JPEG /GrayACSImageDict << /QFactor 0.15 /HSamples [1 1 1 1] /VSamples [1 1 1 1] >> /GrayImageDict << /QFactor 0.15 /HSamples [1 1 1 1] /VSamples [1 1 1 1] >> /JPEG2000GrayACSImageDict << /TileWidth 256 /TileHeight 256 /Quality 30 >> /JPEG2000GrayImageDict << /TileWidth 256 /TileHeight 256 /Quality 30 >> /AntiAliasMonoImages false /CropMonoImages true /MonoImageMinResolution 1200 /MonoImageMinResolutionPolicy /OK /DownsampleMonoImages true /MonoImageDownsampleType /Bicubic /MonoImageResolution 1200 /MonoImageDepth -1 /MonoImageDownsampleThreshold 1.08250 /EncodeMonoImages true /MonoImageFilter /CCITTFaxEncode /MonoImageDict << /K -1 >> /AllowPSXObjects false /CheckCompliance [ /None ] /PDFX1aCheck false /PDFX3Check false /PDFXCompliantPDFOnly false /PDFXNoTrimBoxError true /PDFXTrimBoxToMediaBoxOffset [ 0.00000 0.00000 0.00000 0.00000 ] /PDFXSetBleedBoxToMediaBox true /PDFXBleedBoxToTrimBoxOffset [ 0.00000 0.00000 0.00000 0.00000 ] /PDFXOutputIntentProfile (None) /PDFXOutputConditionIdentifier () /PDFXOutputCondition () /PDFXRegistryName (http://www.color.org) /PDFXTrapped /Unknown /CreateJDFFile false /SyntheticBoldness 1.000000 /Description << /ENU (Use these settings to create PDF documents with higher image resolution for high quality pre-press printing. The PDF documents can be opened with Acrobat and Reader 5.0 and later. These settings require font embedding.) /JPN /FRA /DEU /PTB /DAN /NLD /ESP /SUO /NOR /SVE /KOR /CHS /CHT /ITA >> >> setdistillerparams << /HWResolution [2400 2400] /PageSize [595.000 842.000] >> setpagedevice