Aftershocks and Omori's law in a modified Carlson-Langer model with nonlinear visco-elasticity
aa r X i v : . [ phy s i c s . g e o - ph ] M a y Aftershocks and Omori’s law in a modified Carlson-Langer model with nonlinearvisco-elasticity
Hidetsugu Sakaguchi and Kazuki Okamura
Department of Applied Science for Electronics and Materials,Interdisciplinary Graduate School of Engineering Sciences,Kyushu University, Kasuga, Fukuoka 816-8580, Japan
A modified Carlson-Langer model for earthquakes is proposed, which includes nonlinear visco-elasticity. Several aftershocks are generated after the main shock owing to the damping of theadditional visco-elastic force. Both the Gutenberg-Richter law and Omori’s law are reproduced ina numerical simulation of the modified Carlson-Langer model on a critical percolation cluster of asquare lattice.
PACS numbers: 05.45.-a, 91.30.Dk, 62.20.-x, 45.90.+t
I. INTRODUCTION
Earthquakes follow various types of power laws. Gutenberg and Richter found that the slip-size of faults or theseismic moment S obeys a power law: [1] P ( M ) ∼ − bM , (1)where the magnitude M is defined as M = (2 /
3) log S . The exponent b takes a value ranging from 0.8 to 1.2, andthe most typical value is b = 1. After a large main shock, many aftershocks occur. Omori found that the frequencyof aftershocks decays in a power law: [2, 3] n ( τ ) = K ( c + τ ) p , (2)where τ is the elapsed time from the main shock, and the exponent p is expected to take a value between 0.9 and 1.4.In the original Omori’s law, p is set to be 1. Omori’s law implies that the frequency of aftershocks decays slowly oraftershocks can occur even after a rather long time following the main shock. The mechanisms of earthquakes andaftershocks have been studied by many authors, but are not completely understood. There are several deterministicmodels for earthquakes. The power laws of earthquakes might be understood from the chaotic dynamics of thenonlinear equations. Burridge and Knopoff proposed a block-spring model for earthquakes [4]. Later, Carlson andLanger used a type of velocity-weakening friction and found a power law in the small magnitude range [5]. Manyauthors studied intensively the Carlson-Langer model and its modification. Vieira, Vasconcelos, Nagel pointed out atransition in the magnitude distribution by changing a parameter of the velocity-weakening friction [6]. On the otherhand, there are various models in which the power laws originate from the fractal geometry of faults or asperitiesof plate boundaries. The size distribution of earthquakes has been explained by several authors using such fractalgeometry [7]. We studied the block-spring system on critical percolation clusters with fractal geometry, and showedthat the exponent of the power law is about 1, which is close to the exponent in the original Gutenberg-Richter law [8].There are also several stochastic cellular automaton models for earthquakes. For example, Bak and Tang proposed acellular automaton model for earthquakes which exhibits self-organized criticality [9]. Olami, Feder, and Christensenproposed a modified cellular automaton model with some dissipation [10].The difference between main shocks and aftershocks is not clear in the Carlson-Langer model or the Bak-Tangmodel. Ito and Matsuzaki proposed a modified model of Bak-Tang model and reproduced Omori’s law for aftershocksby randomly perturbing certain regions that slipped during the main shock [11]. Dieterich studied time-dependentfriction in rocks and suggested a relationship between the visco-elasticity and Omori’s law [12]. Hainzl et al. includedsome visco-elasticity in the Olami-Feder-Christensen model and reproduced a version of Omori’s law [13]. Jagla etal. studied a depinning model of viscoelastic interfaces and reproduced a type of Omori’s law [14, 15]. However,the physical meaning of models constructed by Ito-Matsuzaki and Hainzl et al. remains unclear. Jagla’s depinningmodel includes randomness and the relationship between the depinning model and earthquakes is not obvious. In thispaper, we study a modified version of the Carlson-Langer model with visco-elastic elements, and discuss Omori’s lawof aftershocks. Our model is based on mechanics or Newton’s equation of motion, and does not include randomness orexternal noise in contrast to the models by Hainzl et al. and Jagla et al. The time evolution is continuous during slipevents, but stick intervals between successive slips are theoretically evaluated and numerical simulation for the stickintervals is skipped. A rather long-time numerical simulation becomes possible owing to the reduction in computationtime. uy v (cid:129) ¤ (cid:129) ¤ f FIG. 1: Schematic figure of a modified Carlson-Langer model with visco-elasticity.
II. THE CARLSON-LANGER MODEL WITH VISCOELASTICITY
The Carlson-Langer model is a block-spring model. In a one-dimensional model, N blocks of mass m = 1 arecoupled with the nearest neighbor blocks using a spring with a spring constant k . In the original model, each block isfurther coupled with a spring to a rigid plate moving with a constant velocity v . The model equation is expressed as d x i dt = k ( x i +1 − x i + x i − ) + ( v t − x i ) − φ ( v i ) , (3)where v i denotes the velocity of the block dx i /dt , k is the spring constant between neighboring blocks, and φ ( v i )represents the kinetic friction. The maximum static friction is set to be 1. In the Carlson-Langer model, the kineticfriction φ ( v ) is assumed to be φ ( v ) = 1 − σ αv/ (1 − σ ) (4)for v >
0. Two parameters α and σ characterize the velocity-weakening friction. The parameter σ represents thedifference between the maximum static friction 1 and the kinetic friction at the limit of v = 0. The parameters k and σ are fixed to be k = 16 and σ = 0 .
01 in our numerical simulations for the sake of simplicity. The parameter α expresses the velocity-weakening rate of the kinetic friction, i.e., dφ ( v ) /dv = − α at the limit of v = 0. When v i decreases to 0, the slip state changes to the stick state. The stick state is assumed to be maintained until the appliedforce reaches the maximum static friction 1. It is assumed that the reverse motion does not occur [5].In the modified model, a dashpot and spring are added in parallel to the spring coupled with the rigid plate pulledby velocity v as shown in Fig. 1. An additional force f i ≤ i th block by the additional dashpotand spring after the block starts to slip. The additional force works as a braking force to the slipping block. Variousmodels including viscoelasticity are possible. For example, the dashpot and the spring can be added in parallel tothe spring connecting neighboring blocks. These models might be interesting but we investigate only the simplestmodel shown in Fig. 1 in this paper. One reason is that the role of the viscoelasticity as a brake force is rather clearand anther reason is that the dynamics in the stick phase is simple and analytically calculated as explained in laterparagraphs.The force satisfies f i = Ky i where − y i represents the contraction of the additional spring. For a linear dashpot,the force f i satisfies f i = − η ( − u i ) where u i < f i satisfies f i = − η ( − u i ) γ . A power-law relation is assumed between the slip velocity and resistance force, and γ = 1 correspondsto the linear dashpot. A power-law model for the relationship between shear velocity and shear stress is one of thetypical models for non-Newton fluids. The power-law model of 0 < γ < u i is expressed as u i = − ( − f i /η ) β , where β = 1 /γ . Thevelocity difference between the upper plate and the i th block is written as v − dx i /dt which is equal to dy i /dt + u i .Then, x i and f i satisfy d x i dt = k ( x i +1 − x i + x i − ) + ( v t − x i ) + f i − φ ( v i ) ,df i dt = K dy i dt = K ( v − dx i dt − u i ) = K (cid:18) v − dx i dt (cid:19) + K (cid:18) − f i η (cid:19) β . (5)We further assume that v is infinitesimally small, and v is set to be zero during the slip phase. In the stationarystick phase, f i = 0, and v i = 0. First, we calculate the maximum value F of F i = k ( x i +1 − x i + x i − ) + ( v t − x i )among the N blocks. If the i th block takes the maximum value, the time is advanced from t to t + (1 − F ) /v . Thepulling force at the i th block exceeds the maximum static friction 1, and the i th block starts to slip. The neighboringblocks can slip together owing to the elastic coupling of spring constant k . The time evolution of the slip processobeys Eq. (5). Because we assume v is infinitesimally small, two slip events do not occur simultaneously in distantplaces. When the first slip event (main shock) starts, dx i /dt > i , and f i becomes negative. Thatis, the braking force f i works for the slipping block. After the main shock ends and the stick phase starts, the slipvelocity dx i /dt becomes zero, then f i < K ( − f i /η ) β > f i ’s decay to zero. The periodfrom the start of a main shock and the end of the sequence of aftershocks is one shock event. At the end of the shockevent, F i = k ( x i +1 − x i + x i − ) + ( v t − x i ) is smaller than 1 and f i = 0 for all i ’s. Then, we again calculate themaximum value F of F i = k ( x i +1 − x i + x i − ) + ( v t − x i ) among the N blocks, and the upper plate is slowly moveduntil another block begins to slip. The time interval ∆ t until the starting time of the next main shock is calculatedby v ∆ t = 1 − F . Similarly, we can calculate the start time of aftershocks in the sequence of aftershocks. In the stickphase after an aftershock, f i ( t ) can be explicitly expressed as {− f i ( t ) } − β = {− f i ( t ) } − β + (1 − β ) K ( t − t ) η β , (6)for β = 1 because v = 0 and dx i /dt = 0 after the final time t of the aftershock. From Eq. (6), a possible start time t i of the next aftershock for the i th block is evaluated as t i = t + η β K (1 − β ) { ( F i − − β − ( − f i ( t )) − β } , (7)because F i + f i ( t ) = 1 at the start time t i of the next aftershock. The minimum value t s of t i is the actual start timeof the next aftershock, and f i ( t s ) at t = t s is evaluated using Eq. (6). If all F i ’s are smaller than 1, the sequence ofaftershocks ends. In the case of the linear visco-elasticity: γ = β = 1, Eqs.(6) and (7) are replaced by f i ( t ) = f i ( t ) e − K ( t − t ) /η , (8)and t i = t − ηK ln F i − − f i ( t )) . (9)The equation of motion Eq. (5) with Eq. (4) is used only in the slip phases, and the time intervals in the stick phasesbetween main shocks or between two aftershocks are theoretically calculated using Eq. (7) or Eq. (9). We performed anumerical simulation repeating these processes. A characteristic of our model is that a sequence of aftershocks occursafter each main shock and the difference between main shocks and aftershocks is clear. Because v is assumed tobe infinitesimally small, the time of a sequence of aftershocks is infinitesimally small compared with the interval ofsuccessive main shocks, which is an order of 1 /v . III. NUMERICAL RESULTS IN ONE DIMENSION
First, we show the numerical results of the modified Carlson-Langer model of ten blocks to understand the roleof aftershocks. Figure 2(a) shows the time evolution of the sum S = P s i of slip distances at γ = 1 / β = 2), α = 1 . , K = 4, and δ = K/η β = 1. In the time scale of O (1 /v ), aftershocks occur at the same time as the mainshock. The vertical lines show aftershocks. The largest value of S at each slip time corresponds to the main shock.A slip event without aftershocks is shown as a rhombus without vertical lines. The number of rhombi at each sliptime is the sum of one (main shock) plus the number of aftershocks. The number of aftershocks takes a various valueranging from 0 to 6. Figure 2(b) shows time evolution of the profile of x i for the ten-block system using the sameparameters. The light solid lines denote the profiles of x i after the main shocks and the dark dashed lines denotethe profiles of x i after the aftershocks. Figure 2(c) is the enlargement of Fig. 2(b) for 424 . < x < . x i +1 − x i + x i − is positive and large. As the braking forces f i weaken, aftershocks tend to occur near these dentstructures. Aftershocks make the profile of x i flat. When the profile of x i becomes sufficiently flat, a simultaneous x i0.00010.0010.010.11424.25 424.275 424.3 S v t (a) (b) (cid:13)(cid:13) x (cid:13) i(cid:13) (c) FIG. 2: (color online) (a) Time evolution of the sum S of slips for ten-block systems at γ = 1 / β = 2), α = 1 . , K = 4, and δ = K/η β = 1. (b) Time evolution of the profile of x i for ten-block systems using the same parameters. The light solid linesdenote the profiles of x i after main shocks and the dark dashed lines denote the profiles of x i after aftershocks. (c) Enlargementof (b) for 424 . < x < . (a)(c) (b)(d) (cid:129)[ (cid:129)[6 (cid:131)(cid:209)(cid:131)(cid:209)(cid:131)(cid:209) (cid:131)(cid:209) nn nn FIG. 3: (a) Semi-logarithmic plot of the probability distribution of the elapsed time τ of aftershocks at γ = 1 , α = 1 , K = 0 .
5, and δ = K/η = 0 .
05. (b) Double-logarithmic plot of the probability distribution of τ at γ = 1 / , α = 1 , K = 4, and δ = K/η = 1.(c) Double-logarithmic plot of the probability distribution of τ at γ = 1 / , α = 0 . , K = 0 .
4, and δ = K/η = 10000. (d)Double-logarithmic plot of the probability distribution of τ at γ = 1 / , α = 0 . , K = 0 .
4, and δ = K/η = 100000. slip or a large shock occurs next. Through chaotic dynamics in the simultaneous slip event, the spatial fluctuationincreases over time and another dent structure appears after the large-scale shock. Similar intermittent dynamicswere also observed in the original Carlson-Langer model [8]. In the original Carlson-Langer model, small slip evens(small-scale main shocks) frequently occurs near the dent structures. In our modified model, aftershocks play a rolein these small slip events in the original Carlson Langer model.We investigated the statistical properties of aftershocks in a one-dimensional system of N = 1000. We performed P (cid:13) M(cid:13) 0.001(cid:13)0.01(cid:13)0.1(cid:13)1(cid:13)-2(cid:13) -1(cid:13) 0(cid:13) 1(cid:13) P (cid:13) M(cid:13) 0.001(cid:13)0.01(cid:13)0.1(cid:13)1(cid:13)-2(cid:13) -1(cid:13) 0(cid:13) 1(cid:13) P (cid:13) M(cid:13) (a) (b) (c)
FIG. 4: Probability distribution of the magnitude at (a) α = 0 .
3, (b) 0.4, and (c) 0.5 for γ = 1 / , K = 0 .
4, and δ = K/η =100000. The dashed line in Fig.4(b) denotes a curve of P ( M ) ∝ − . M . numerical simulations by changing parameters γ, α, K and δ = K/η /γ . We present typical examples in this paper.Figure 3(a) shows a semi-logarithmic plot of the probability distribution of the time τ when aftershocks occur after themain shocks at γ = 1 , α = 1 , K = 0 .
5, and δ = K/η = 0 .
05. The dashed line is n ( τ ) = 0 .
048 exp( − . τ ), althoughthe dashed line overlaps with the solid curve and is hardly visible. In this case of linear visco-elasticity γ = 1,the probability distribution obeys an exponential distribution. This implies that the simple linear model does notreproduce Omori’s law of aftershocks. Figure 3(b) shows a double-logarithmic plot of the probability distribution of τ at γ = 1 / , α = 1 , K = 4 , and δ = K/η = 1. The dashed line is n ( τ ) = 356 / (373 + τ ) . Figure 3(c) shows a double-logarithmic plot of the probability distribution of τ at γ = 1 / , α = 0 . , K = 0 . δ = K/η = 10000. The dashedline is n ( τ ) = 0 . / (31 + τ ) . . Figure 3(d) shows a double-logarithmic plot of n ( τ ) at γ = 1 / , α = 0 . , K = 0 . δ = K/η = 100000. The dashed line is n ( τ ) = 0 . / (310 + τ ) . . Similarly, the probability distribution at γ = 3 / , α = 1 , K = 1 , δ = K/η /γ = 0 . n ( τ ) = 60000 / (350 + τ ) . is obtained. The probability distributions approximately obey Omori’s law in the nonlinear visco-elastic models. Theexponents p are evaluated at p = 2 . , , . .
02 for γ = 3 / , / , / /
4, respectively. The exponent p isinterpreted as p = ∞ for the linear case of γ = 1. These results suggest that the exponent p decreases with γ . For γ = 1 /
4, the exponent p is close to the exponent p = 1 of original Omori’s law.Figures 4 shows the probability distribution P ( M ) of the magnitude at (a) α = 0 .
3, (b) 0.4, and (c) 0.5 for γ = 1 / , K = 0 .
4, and δ = K/η = 100000. Figure 4(c) is a magnitude distribution for the numerical simulationshown in Fig. 3(d). In the original Carlson-Langer model, the magnitude distribution has a hump structure arounda rather large value of M at α > M at α <
1. The magnitude distibution exhibits apower law behavior near α = 1. Even in our model with nonlinear viscoelastisity, the magnitude distribution decaysrapidly for large M at small α as shown in Fig. 4(a) and has a long tail at large α as shown in Fig.4(c), althoughthere is a shoulderlike structure near M = − . α = 0 . M > − .
5, although the linearity is not so good as compared to the modelson a percolation cluster as shown in Fig. 6(a) and Fig. 7(a). The exponent b is evaluated as 1.2, which is close to the b value of the Gutenberg-Richter law for earthquakes. IV. NUMERICAL RESULTS ON A PERCOLATION CLUSTER
Next, we investigated the statistical properties of aftershocks on a critical percolation cluster on a square latticeof 150 × k = 16. We assumed that blocks are locatedonly in the largest mutually-connected percolation cluster. In the previous paper, we performed a numerical simulationof the original Carlson-Langer model on the same critical percolation cluster [8]. We show numerical results for themodified Carlson-Langer model on the same percolation cluster at γ = 1 / α = 1 . , K = 1, and δ = K/η = 0 . S of the main shocks. Figure 5(b) shows the time evolutionof the aftershocks after a main shock at v t = 3 . S (a) N S (b)(c) (d) (cid:13) v t (cid:131)(cid:209) S j i a FIG. 5: (color online) (a) Time evolution of main shocks at γ = 1 / α = 1 . , K = 1, and δ = K/η = 0 .
05. (b) Time evolutionof aftershocks after the main shock at v t = 3 . S of themain shock and the number N a of aftershocks after each main shock. The dashed line is N a = 9 S . (a) n P (b) (c) (cid:131)(cid:209) (cid:131)(cid:209) n -5 -6 -5 n (cid:13) -6 M<0 M>0
FIG. 6: (a) Probability distribution of the magnitude M = (2 /
3) log S at γ = 1 / α = 1 . , K = 1, and δ = K/η = 0 . P ( M ) ∝ − . M . (b) Probability distribution of τ . The dashed line is n ( τ ) = 13 . / (680+ τ ) . . (c) Probability distribution of τ for M >
M <
0. The dashed lines are n ( τ ) = 7 . / (2560+ τ ) . and n ( τ ) = 598 / (916 + x ) . . for each slip event. Figure 5(d) shows a relationship between the total slip size S of the main shock and the number N a of aftershocks after each main shock. The number of aftershocks increases with the magnitude of the main shock.The dashed line is N a = 9 S . .Figure 6(a) is the probability distribution of the magnitude M = (2 /
3) log S at γ = 1 / α = 1 . , K = 1, and δ = K/η = 0 .
05. The dashed line denotes a distribution P ( M ) ∝ − . M . A good agreement is seen for M < .
5. Atype of Gutenberg-Richter law is satisfied. Figure 6(b) shows the probability distribution of the elapsed time τ of the (a) (b) (cid:131)(cid:209) n P -5 n (cid:13) (cid:131)(cid:209) -5 -6 (c) M>0M<0
FIG. 7: (a) Probability distribution of the magnitude M = (2 /
3) log S at γ = 1 / α = 0 . , K = 0 .
1, and δ = K/η = 100000.The dashed line denotes an exponential distribution P ( M ) ∝ − . M . (b) Probability distribution of τ . The dashed line is n ( τ ) = 0 . / (48+ τ ) . . (c) Probability distribution of τ for M >
M <
0. The dashed lines are n ( τ ) = 0 . / (112+ τ ) . and n ( τ ) = 0 . / (215 + x ) . . aftershocks after the main shocks. The dashed line shows n ( τ ) = 13 . / (680 + τ ) . . Omori’s law of p = 1 . τ after the main shocks of larger magnitude M >
M <
0. Figure 6(c) shows the two probability distributions. Both distributions obey Omori’slaw but the exponent p is 1.35 for M >
M <
0. We also performed a numerical simulation on the samepercolation cluster on a square lattice of 150 ×
150 at γ = 1 / , α = 0 . , K = 0 .
1, and δ = K/η = 100000. Figure7(a) is the probability distribution of the magnitude M = (2 /
3) log S . The dashed line denotes an exponentialdistribution P ( M ) ∝ − . M . A good agreement is observed for M < .
5. The distribution P ( M ) ∝ − bM with b = 1 corresponds to the original Gutenberg-Richter law. Figure 7(b) shows the probability distribution of τ . Thedashed line shows n ( τ ) = 0 . / (48 + τ ) . . Omori’s law of p = 1 .
02 is obtained. The Gutenberg-Richter law andOmori’s law close to the original versions are obtained at this parameter set. Figure 7(c) shows the two probabilitydistributions of τ for M >
M <
0. Both distributions obey Omori’s law but the exponent p is 1.1 for M >
M <
0. We do not understand the reason but the exponent seems to depend on the magnitude of the mainshock and the exponent p is larger for larger M in our model. We also performed a similar numerical simulation for alinear visco-elastic model of γ = 1 at α = 1 , K = 0 .
5, and δ = K/η = 1. We found that the probability distribution of τ obeys an exponential function in the linear model, which is the same as the behavior of the one-dimensional modelshown in Fig. 3(a). These numerical results suggest that the exponent p of Omori’s law does not strongly depend onthe dimension. V. SUMMARY
We have proposed a modified Carlson-Langer model which can generate aftershocks, and performed numericalsimulations. The visco-elastic element works as a brake in the slip process of blocks. Owing to the damping ofthe braking force, the forward tensile force for blocks increases, and aftershocks are induced. By using nonlinearvisco-elastic elements, Omori’s law has been reproduced. On the other hand, Omori’s law could not be reproducedin simple linear visco-elastic models. Both the Gutenberg-Richter law with exponent close to 1 and Omori’s law withexponent close to 1 are observed in the numerical simulation of the modified Carlson-Langer model with γ = 1 / [1] B. Gutenberg and C. F. Richter, Seismicity of the Earth and Related Phenomena (Princeton University Press, Princeton,NJ, 1954).[2] F. Omori, J. Col. Sci. Imp. Univ. Tokyo , 111 (1894).[3] T. Utsu, Geophysical Magazine , 521 (1961).[4] R. Burridge and L. Knopoff, Bull. Seismol. Soc. Am. , 341 (1967).[5] J. M. Carlson and J. S. Langer, Phys. Rev. Lett. , 2632 (1989). [6] M. de Sousa Vieira, G. L. Vasconcelos, and S. R. Nagel, Phys. Rev. E , R2221 (1993).[7] V. De Rubeis, R. Hallgass, V. Loreto, G. Paladin, L. Pietronero, and P. Tosi, Phys. Rev. Lett. , 2599 (1996).[8] H. Sakaguchi and S. Morita, J. Phys. Soc. Jpn. , 114006 (2013).[9] P. Bak and C. Tang, J. Geophys. Res. , 15635 (1989).[10] Z. Olami, H. J. S. Feder, and K. Christensen, Phys. Rev. Lett. , 1244 (1992).[11] K. Ito and M. Matsuzaki, J. Geophys. Res. , 6835 (1990).[12] J. H. Dieterich, J. Geophys. Res. , 3690 (1972).[13] S. Hainzl, G. Z¨oller, and J. Kurth, J. Geophys. Res. , 7243 (1999).[14] E. A. Jagla, F. P. Landes, and A. Rosso, Phys. Rev. Lett. , 174301 (2014).[15] E. A. Jagla, Phys. Rev. E90