The mechanisms of spatial and temporal earthquake clustering E. A. Jagla1 and A. B. Kolton1 1Centro At´omico Bariloche, Comisi´on Nacional de Energ´ıa At´omica, 8400 Bariloche, Argentina∗ The number of earthquakes as a function of magnitude decays as a power law. This trend is usually justified using spring-block models, where slips with the appropriate global statistics have been numerically observed. However, prominent spatial and temporal clustering features of earthquakesarenotreproducedbythiskindofmodeling. Weshowthatwhenaspring-blockmodel iscomplemented with amechanism allowing for structuralrelaxation, realistic earthquakepatterns 9 areobtained. Theproposedmodeldoesnotneedtoincludeaphenomenological velocityweakening 0 friction law, as traditional spring-block models do, since this behavior is effectively induced by the 0 relaxational mechanism as well. In thisway,themodel providesalso a simple microscopic basis for 2 thewidely used phenomenological rate-and-state equations of rock friction. n a I. INTRODUCTION responsible for less than about 5 % of the total released J seismic moment, the finding that the GR law is also 3 obeyed within individual aftershock sequences strongly 1 The distribution of earthquakesin nature follows non- suggeststhat consistentandcompatible explanationsfor trivial patterns, some of which are captured by well ] GR and Omorilaws should exist. We may thus say that n knownempiricallaws. The Gutenberg-Richter (GR) law at present, there is not a single, unified picture of the n [1,2]statesthatthenumberofearthquakesasafunction physics behind some of the most robust features of seis- - of magnitude N(M) scales as N(M) ∼ 10−bM. The ex- s micity, namely GR andOmorilaws. In addition, the use i ponentbisverynearly1. TheOmorilawreferstotempo- d of rate-and-state equations, although widely supported ralcorrelationsbetween earthquakes,in particular to af- . by experimentalresults,remainsessentiallya crude phe- t tershocks,namelythetemporalclusteringofearthquakes a nomenologicalapproach. m followingalargeone,usuallycalledthe mainshock. The Weshowherethatwhenaspring-blockmodelwithout Omorilawofaftershocks[7]statesthatthenumberofaf- d- tershocksper unitof time decaysas∼(t+c)−p with the any a priori velocity dependent frictional force is com- plemented with an appropriate relaxational term as dis- n time t from the main shock. The exponent p is typically cussedbelow,itproduces: 1)earthquakepatternsandin o very close to 1, and c is a time constant of the order be- c particular aftershock sequences quantitatively compara- tween minutes and hours. Aftershocks occur mainly in [ ble with real ones, 2) a velocity weakening friction law, the spatial region where the rupture of the main shock and in general, agreement with the predictions of the 1 took place. rate-and-state equations, and 3) a power law decay of v 7 The GR law has been shown to be compatible with number of earthquakes with magnitude compatible with 0 a state of (at least partial) critical organization of the the GRlaw,withanexponentb thatcompareswellwith 9 system[3, 4, 5]thatis understandableinterms ofspring- actual values. 1 block models, when an appropriate velocity weakening 1. friction law (i.e., a friction force decreasing with the 0 relative velocity of the sliding elements) is assumed to II. MODEL IN THE ABSENCE OF 9 hold. This kind of modeling was pioneered by Bur- 0 ridge and Knopoff[6] (BK), and was extended along dif- RELAXATION : v ferent directions afterwards, particularly in the works i on self-organized-criticality of the eighties [3, 4]. The Our modeling is based on the original BK model[6], X BK model reproduces the global statistical behavior im- with the important difference that the friction law be- r plied by the GR law, but it fails to account for the tween the blocks and the substrate is not a priori as- a existence of spatial and temporal correlations observed sumed to have any particular form such as the velocity- in actual seismicity. On searching for the origin of the weakening form commonly used. A velocity dependent aftershock phenomenon, Dieterich[8] followed by others friction law emerges naturally at large scales from the [9,10,11,12]haveshownthatananalysisbasedonrate- characteristic collective dynamics of elastic manifolds in and-state equations[13, 14] is able to justify the appear- random media[16]. In this context, an elastic interface ance of aftershocks following the Omori decay. On this (whichcorrespondstotheblocksjoinedbyspringsinthe perspective, it is puzzling that the use of a BK model BKmodel-wealreadydescribethetwodimensionalcase, with a rate-and-state friction law does not produce re- more appropriate to real faults) is driven through a dis- alistic aftershocks[15]. Although aftershocks usually are orderedpotentialenergylandscape(the‘substrate’)that models the random nature of asperities. The velocity- dependent frictional force between the blocks and the substrateofthe BKmodelis thereforereplacedbya dis- ∗E-mail: [email protected] ordered potential energy landscape that is chosen ran- 2 domly and uncorrelated for each point block of the dis- (A) cretizedinterface. Inconcrete,ourmodelisdescribedby the overdamped equation ∂u λ i,j =k0∇2ui,j +fi,j +k1(X0(t)−ui,j) (1) ∂t where u is a continuous variable representing the dis- i,j placement of the block labeled by the indices (i,j) in a two dimensional grid, X0(t) is the driving variable (usually X0(t) = Vt, we will refer also to X0 and to k1(X0(t)−ui,j) as the strain and local stress in the sys- (B) (C) tem,respectively),∇2 isthediscreteLaplacianoperator, and f is the pinning force at each block, which is as- i,j sumedtobeshort-rangecorrelatedalongthedirectionof aftershock the block displacements. For numerical convenience this epicenter main shock spatially random force is chosen in the following way: a epicenter random position u0 is selected for each i,j, then the i,j force is f =k(u −u0 ). When, upon the dynamical i,j i,j i,j evolution, f reaches some threshold value fth (that is i,j i,j chosenrandomlydistributedbetween1+κand1−κ),the corresponding u0 is given the value u +δ, where δ is i,j i,j randomlychosenbetween −1and 1. Also,a new thresh- FIG.1: Onedimensionalsketchofmainprocesses thatoccur oldisassignedtositei,j. Similarresultsareobtainedby in our model. (A) In the absence of relaxation, the inter- using standard, though computationally more demand- face(verticalline,withcoordinatesui,j)isdriventotheright ing, short-ranged correlated smooth pinning forces such by an external force (not shown) through a set of randomly as the ones used in Refs. [17, 18]. placed pinning centers, represented by the gray rectangles. The mid point of the rectangles have coordinates u0 , and i,j We use periodic boundary conditions, and from now the horizontal length is the range in which pinning is effec- on, we set the values k0 = 0.1, k = 1, κ = 0.8. Taking tive. Differentlengthsindicatedifferentthresholdvaluesfth. i,j also the numerical lattice constant as unity, this renders Upon driving, the system passes from the configuration in- our time, distance, and forces, dimensionless. We have dicated by the continuous line to the dashed line, and an tried other parameter sets, finding no qualitatively new event is triggered at the site indicated by the arrow, where results. ThevalueofλinEq. (1)fixesthetimescalenec- the maximum local pinning force is overpassed. The system essaryforthesurfacetoadapttothe conditionsdictated goes through a cascade process (not fully indicated) onto a by u0 and X0. We work in the case λ → 0, i.e., we give new meta-stable configuration (dotted line) in which some pinningcenters havebeen refreshed (outlined rectangles). In the surface time to relax to what we call a meta-stable (B)structuralrelaxationisacting. Aquiterelaxed(andthere- configuration,defined by equating the right hand side of foremorecoherentlypinned)initialconfiguration(continuous Eq. (1) to zero, for fixed values of X0 and u0i,j. Note line) is driven until a main shock occurs, at the configura- in particular that this also means that the duration of tion corresponding to the dashed line. After the system has individual earthquakes is zero in our implementation. reachedameta-stableconfiguration(dottedlineandoutlined rectangles) relaxation continuous to act modifying the posi- Abrupt rearrangements occur in the system whenever tion of the pinning centers. The arrows at the center of the atsomeparticularpositioni,j theforcefromthepinning rectanglesindicatethelocalvalueof(u−u0). Thearrowsjust potentialofthe substrateisnotabletosustainanymore outsidetherectanglesindicatethevaluesofdu0/dtaccording the surface pinned to it (see an example of this situation toEq. (2),thatproduceadriftinthepositionofthepinning in Fig. 1A). In this situation the local rearrangement centers. This drift may cause a further instability as can be of the surface can trigger instability events in neighbor seen in (C). The process in (C) is an aftershock to the main sites, and the process continues until the surface finds a shock in (B). newgloballystableconfiguration. Thefullsequenceofall rearrangementstriggeredby aninitialinstability iswhat we call an event, or an earthquake,being the position of A. Results the triggering instability its epicenter. We measure the seismic moment m0 of events as m0 = Pi,j∆i,j where ∆ is the displacement caused by the event at position Ourmodelintheabsenceofrelaxationhasbeenwidely i,j i,j. Inordertocomparewithrealearthquakes,themag- studied, and is known to have a well defined size dis- nitudeM ofaneventisdefinedfromtheseismicmoment tribution of events[17, 18, 19, 20] (note that the usual as M =2/3log10m0. definition of the decay exponent τ is given as a func- 3 FIG. 2: Results without relaxation. (A) Magnitude his- togram for systems of 512x512 sites, with k1 = 0.01 (full symbols) and k1 = 0.001 (open symbols). Continuous line hasaslope b=0.4. (B)Magnitude-time plot fora systemof FIG. 3: Results for the model with relaxation as compared size 256x256, k1 =0.01. toearthquakesinCalifornia[24]. (A)Magnitude-timeplotfor a 512x512 system in the presence of relaxation (k1 = 0.01, R/V =500). (B) Actualearthquakesin theCalifornia area. tion of our b as τ = 2b/3+ 1). In Fig. 2A we show this distribution. Weobserveapowerlawwithexponent b ≃ 0.4, which is consistent with the expected results sion of a simple additional ingredient changes this sce- for two dimensional elastic interfaces both from scaling nario drastically. This ingredient turns out to be what arguments[19] and from recent analytical and numerical we have called structural relaxation [21]. The primary calculations[20]. Thisvalueishoweverwelldifferentfrom the b≃1 observed for actual earthquakes. An exponen- physical justification of its inclusion is the following. It is known that in solid friction the friction coefficient at tialcut-off for largeeventsize exists due to confinement. rest increases with the time the surfaces have been in This cut-offiscontrolledbythe rigidityk1 ofthe driving contact [22, 23]. This is telling us about the existence of spring and occurs when the spatial extent of the events in the direction of the displacements is of order k−ζ/2, a temporaldependent mechanismthat makesthe sliding 1 surfacesgetmore attachedorpinned toeachotherwhen withζ theinterfaceroughnessexponentatlowvelocities. they remain in contact for a longer period of time. In The crossoverto the exponentialbehaviorthus moves to the present model, a rather simple way to include such larger magnitudes as k1 is decreased [17, 18, 19, 20]. aneffectistoconsider,inadditiontotherandompinning InFig. 2Bweshowatime-magnitudeplotofallevents force that the substrate performs on the sliding surface, occurringinaparticulartimeinterval. Itisapparentthat thereactionthatthesurfaceperformsontothesubstrate. no obvious temporal correlations occur in this case, and If we give the substrate the possibility to react to this morequantitativeobservationsconfirmthis fact. Spatial force,the systemwillgainpinning energyby makingthe correlations are not observed neither. This is qualita- substratemorecorrelated,sotopinbetterthecorrelated tively similar to what has been obtained for the original interface structure. This is in general a slow process, BK model[6]. and the longer the surface remains in contact with the substrate, the stronger the join. This process of attach- mentis howeverstoppedandrestartedwhena slipevent III. MODEL IN THE PRESENCE OF occurs, since the values of the disordered potential re- RELAXATION freshen, becoming uncorrelated again. We will show in thefollowingthatthissimplemechanismisenoughtoex- Sofar,inthepresentform,themodeldoesnotgiveany plain, in particular,the appearance of a robust sequence clueonthereasonforearthquakeclustering,ortheorigin of aftershocks, and the occurrence of a velocity weaken- of a velocity weakening friction law. However, the inclu- ing friction law. Our structural relaxation mechanism is 4 D t+c [ days] 2 0.1 1 10 100 1000 104 N(t) 103 102 101 100 1E-4 1E-3 0.01 0.1 1 D t+c 1 [t] FIG.4: Spatialdistributionofaftershocksinthesimulations. The slip surface of a large event (shadowing proportional to FIG. 5: The Omori law. Histogram of the number of events the local slip) and the before- (full) and after-events (open aftermainshocksinthesimulation(opensymbols),averaged symbols) of magnitude larger than 0 occurring in a symmet- over 120 main shocks, and the histogram of aftershocks for rictimeintervalδt=0.0022aroundthemaineventareshown. events in the California area (full symbols), averaged over 7 The increase in seismicity after the main shock is clearly ob- events of magnitude M > 6.0 in the time period considered. servable, as well as the localization of aftershocks mainly at ∆tis thetimesince main shock. Curveshavebeen vertically andneartheslipsurfaceofthemainshock. Thefiguredepicts rescaled, setting the value 1 for large ∆t. Continuous lines a portion of size 350x350 of a system of 512x512. are fittings to Omori law with p = 1. The time shifts are c1 =3×10−4, c2 =0.05 days. what in other contexts is called the ageing of the mate- rial. afterthemainshock,maydestabilizeandgenerateanew The modification to the model is as follows. We allow event (Fig. 1(C)). Note that in this description it is ob- the values of u0 to relax in time according to i,j vious that aftershockswill occur at, or near, the rupture regionofthemainshock. Itisalsoworthnotingherethat ∂u0 i,j =−R∇2f =Rk∇2(u0 −u ). (2) aftershocks are triggered by the inner dynamics of the ∂t i,j i,j i,j system, and that most aftershocks occur also if we stop thedrivingofthesystemafterthemainshock. Theseem- This conserveddynamicsfor the shiftofthe disorderpo- ingly contradictoryfact that aftershocks (which must be tential at different block points is a generic way to in- triggered by an initial instability) are originated in a re- troduce the back effect of the surface on the substrate (the actualization of the u0’s when a slip event occurs is laxationmechanismisunderstoodwhenonerealizesthat due to the disorder,relaxationaccordingto Eq. (2) may made as before). The coefficient R is a measure of the produce local increases in the forces f , and this can intensityorrateofrelaxation,andcanthusberelatedto i,j triggeraftershocksif a localthresholdis overcome. Note experimental relaxation times. Equation (2) generates in this respect the opposite direction between (u−u0) a tendency for the local forces fi,j to become uniform and du0/dt at the aftershock epicenter, in Fig. 1C. acrossthesystem,generatingastrongercontactbetween surface and substrate. This relaxational effect competes with the driving, which forces the movement of the sur- face onto the substrate at a fixed average velocity V. A. Results and comparison with an actual earthquake sequence The relevant parameter that measures the competition between the two effects is the ratio R/V. Themechanismbywhichearthquakeclusteringoccurs InFig. 3Aweshowamagnitude-timeplotofeventsin canbesummarizedasfollows(seeFig. 1). Ifaparticular asimulationofasystemof512x512sites,inthepresence regionofthe samplehasnotexperiencedalargeeventin of relaxation (R/V = 500). For comparison, the same aratherlongperiodoftime,thestructuralrelaxationhas plot for the earthquakes in the California area [24] is madethisregionstronger(Fig. 1(B)).Whenaneventoc- presented as Fig. 3B. The visual similarity is striking. curs (driven by the overalldisplacementbetween surface Inbothgraphs,the verylargeactivityimmediately after andsubstrate)thecontactsrefreshenandlargevariations largeevents,i.e.,theexistenceofaftershocks,isapparent. in the local forces remain. Relaxation continues to act, It is worth noting here that a minimum value around trying to uniformize the local forces. In this process, R/V ∼100of relaxationis necessaryin order to observe particularpoints that wereoriginallystable immediately aftershocks in the model, so the structural relaxation is 5 nts cummulative number of events (arb.units) A model (arb.units) cummulative seismic moment cummulative number of eve 134650500000200000.93A2 2.934modelt 2.93116.. 01 (arb.units)cummulative seismic momen cummulative number of events (arb.units)2.90B e a r athctquuaal kes t 2.9 (arb.units)5 cummulative seismic moment cummulative number of events255000000 B 6630 e 6a 6 r4 a0tth c(tqduuaaay6l ks65)e0s 899... 901 (arb.units)tcummulative seismic moment 0 3000 6000 9000 t (days) FIG.6: Cumulativenumberofeventsandcumulativeseismic FIG. 7: Detail to Fig. 6. Cumulative number of events and moment for thesequences presented in Fig. 3. cumulative seismic moment (taken as zero just before the main shock), following the events indicated by vertical ar- rows in Fig. 6. In both cases, dotted lines are fitting to the the crucial ingredient behind these particular effects. (cumulative) Omori law with p=1. The spatial location of aftershocks are strongly corre- latedwiththeslipsurfaceofthemainevent. InFig. 4we Theevolutionoftheaccumulatedseismicmomentisalso show the regionthat has slip in a large event in the sim- qualitatively similar in both cases, with the main shock ulations,togetherwiththeepicentersofalleventsoccur- accounting for most of the released seismic moment of ringinasymmetrictimeintervalaroundthemainshock. the whole sequence. Events before and after the main shock are shown in a separate way. We see that there is a rather uniform spa- Finally,inFig. 8wepresentananalysisofthetimein- tial distribution of before-events,whereas once the main tervals∆t betweensuccessiveeventsofmagnitude larger shock has occurred, aftershocks occur at and near the than some defined threshold M0. The main characteris- region in which the main slip occurred. tics observed for the real sequence, that are reproduced by our model are the following. The curves are roughly In Fig. 5 we plot the histogram with the number of events after a main shock as a function of time. In the independent of the threshold value M0 chosen, and the global behavior represents almost an exponential decay simulations, we average over 120 large events in a single with ∆t, but with a reproducible excess of events at low simulation of size 512x512. The continuous lines corre- spondtoOmorilawsoftheformN(t)=A/(t−t0+c)p+ ∆t. This excess is accounted for by the aftershocks. N0, where t0 is the time of the main shockandN0 is the value of background seismicity. For reference, we also plot the number of aftershocks of the seven earthquake IV. AVERAGED FRICTIONAL PROPERTIES with M > 6.0 in the California area in the considered time period. Both cases are well fitted by an Omori law The second set of results that we present corresponds with p∼1. tothestress-strainrelationofthemodelfordifferentdriv- Figure 6 shows plots of cumulated number of events ing protocols. First of all we recall that the model with- andseismicmomentcorrespondingto the sequencespre- out relaxation shows a stress σ that is independent of sented in Fig.3, and Fig. 7 is a detail after the events the strainrate,sincethe internaltime scaleofthe model indicated by vertical arrows in Fig. 6. The cumulative is very rapid compared to the driving. In Eq. (1) this number of events is fitted in both cases with a cumula- means that λ/V → 0. The inclusion of relaxation in- tive Omori law, with p = 1 with very good agreement. troduces a new time scale (set by the parameter R in 6 1.18 s (t) s A 1.17 av 1.16 s (t) 1.15 t 1.14 1.13 1.12 1.12 t 1.10 0.01 0.1 1 10 V/R100 1.22 s 1.20 B 1.20 s p s 1.15 1.18 p 10-4 10-3 10-2 t 1.16 h 1.14 t h 1.12 2.90 2.91 2.92 2.93 2.94 2.95 time FIG.8: Timeintervalsdistribution. Distributionofthetimes FIG. 9: Global frictional properties of the model. (A) Mean (normalizedbythemeantime∆tm)betweeneventsofmagni- stressinasystemof256x256asafunctionofrelativevelocity. tudelarger than M0 for our data (A) and for earthquakesin Adetailofthetemporaldependenceofstressisgivenfortwo California(B).IndependenceofM0,andaratherexponential points,emphasizingthelargerfluctuationsthatappearwhen distribution withan excessduetoaftershocksfor small ∆tis relative velocity is lower. (B) Time evolution of stress in a observed in both cases. systeminwhichvelocityischangedfromV/R=0inthehold periods(indicatedbyarrows), toV/R=1intherest oftime (resultsshowncorrespondtoanaverageovertenrealizations). Eq. (2)) and now the average stress in the system de- Inset: The value of the stress peak as a function of the hold pends on the ratio R/V. When R/V is small, the effect time. of relaxation is negligible, and the stress will be similar to that in the absence of relaxation. However, if R/V is high enough, relaxation will act by effectively correlat- ing some time interval (the hold time) and then is re- ing the pinning potential in a larger spatial region. The initiated. First of all, a logarithmic decrease of stress size of this region increases with R/V. A spatially more during the hold time is observed. This occurs because correlated pinning potential produces, in turn, a larger the system continues to relax during the hold time and averagestressinthesystem. Weconcludethatthelarger some instability events continue to occur for some time. is R/V the larger is the average stress. In other words, Thisisrelatedtoourpreviousstatementthataftershocks the model will display velocity weakening. In Fig. 9A also occur if driving is stopped after a main shock. De- we show a plot of the stress in the system as a function spite the stress reduction during the hold time, a stress of strain rate where this weakening is clearly observed. peak occurs after re-initiation of sliding. The height of For large strain rates the stress converges to the value this peak increases logarithmically with the hold time. correspondingto no relaxation,whereasthe behaviorfor This peak is a consequence of the more stable configu- very small strain shows a saturation at a larger value. ration that the system reached due to relaxation during The transition between these two values is logarithmic the hold time. The phenomenon is similar to the one andspans abouta factorof 100ofstrainrate. Note that observed in glass forming materials, where it has been the values reported above as necessary to observe after- explained using the same ideas [21]. These results are in shocks (R/V & 100) correspond to the limit of small remarkable agreementwith those obtained in laboratory velocity in this plot. A closer examinationat the instan- measurements [22, 25]. taneous stress-strain relation reveals that the lower the strain rate, the more pronounced the fluctuations in the instantaneous stress. V. GUTENBERG-RICHTER BEHAVIOR Additionalinformationonthefrictionalbehaviorisob- tainedbystudyingthesystemresponsetoabruptchanges The inclusion of relaxation produces also a change in of the strain rate. We show in Fig. 9B in particular, thedecayingexponentboftheGRlaw. InFig. 10wesee the stress on a system in which driving is stopped dur- that a power law decaying is maintained in the presence 7 value b∼1 observed in actual seismicity. N(M)1 0R / V 0 k.011 Thefactthatthebexponenttakesavaluecloseto1in the presence of relaxation, quite independent of the pre- 100 0.02 0.1 100 0.01 cise value of the relaxation parameter and other details 100 0.005 0.01 500 0.01 of the model seems to indicate that relaxation takes the 2000 0.01 system out of its original universality class with b≃0.4, 1E-3 toanewonewithb≃1.0. Coincidenceofthisvaluewith 1E-4 actual ones is another indication that we are capturing essentialfeaturesoftheseismicprocesswiththeinclusion 1E-50 1 2 3 4 of the relaxation mechanism. M VI. CONCLUSIONS FIG. 10: The decaying of number of events with event mag- Summarizing, in the present paper we have presented nitude. The R= 0 case is included again for reference. The a modeling that combines a spring-block type system in thincontinuouslineshowsthecase b=1,followed ratherac- the spirit of the BK model but without a priori velocity curatelybyactualseismicevents. Theresultsinourmodelin the presence of relaxation show a behavior compatible with weakeningfriction,witharathergenericimplementation the actual seismicity, with a cut-off at large event size that ofageingeffectswithintheslidingmaterials. Themotiva- increases upon decreasing the spring constant k1. However, tionforthisapproachwastointroduce,inaspring-block finite size effects seem to be appreciably for the system sizes model, a mechanism that generates (and not merely as- used. The larger decaying rate observed for M < 1 is an sumes) non trivial frictional effects, which can produce spuriouseffectassociatedwitheventscomparableinsizewith realistic temporal and spatial clustering of earthquakes. our numerical mesh. Our model allows to obtain a time sequence of events that globally follow the GR law with a b ≃ 1 expo- nent,andatthe sametime highlynon-trivialspatialand of relaxation, with a b value substantially larger than temporalcorrelationscompatible, inparticular,with the that corresponding to no relaxation. Once a minimum value of relaxation has been over passed (R/V ∼ 20), Omori law. In addition we have shown that frictional properties of the model compare very well with labora- the b decaying exponent is quite insensitive to the pre- tory results. We think this model gives a unified and cise value of relaxation. There seems to be an excess of comprehensive physical picture of all these phenomena. events of large magnitude, before the cut off is reached. The cut off and the peak corresponding to large events are mainly dependent on the value k1 of the spring driv- ing the system. The smaller this value, the larger is the VII. ACKNOWLEDGMENTS cut-off. It is not clear however if this tendency can be extrapolated to very small values of k1. Unfortunately, ThisresearchwasfinanciallysupportedbyConsejoNa- to simulate decreasing values of this spring constant re- cional de Investigaciones Cient´ıficas y T´ecnicas (CON- quires an increase in system size, and we reach rapidly ICET),Argentina. PartialsupportfromgrantsPIP/5596 very time consuming runs. The obtained decaying expo- (CONICET) and PICT 32859/2005 (ANPCyT, Ar- nent in the presence of relaxationis compatible with the gentina) is also acknowledged. [1] C. H. Scholz, The Mechanics of Earthquakes and Fault- Sci Imp Univ Tokio 7:111-200. ing, (Cambridge University Press, Cambridge, England, [8] DieterichJ-H(1994)Aconstitutivelawforrateofearth- 2002). quakeproductionanditsapplicationtoearthquakeclus- [2] GutenbergB,RichterCF(1956) Magnitudeandenergy tering. J GeophysRes 99:2601-2618. of earthquakes. Ann Geophys (C.N.R.S) 9:1. [9] Marcellini A(1997) Physicalmodelofaftershocktempo- [3] Bak P,TangC, Wiesenfeld K (1987) Self-organized crit- ral behavior. Tectonophysics 277:137-146. icality: An explanation of the 1/f noise. Phys Rev Lett [10] Ziv A, Rubin A M (2003) Earthquakes and Acoustic 59:381-384. Emission. J Geophys Res B: Solid Earth 108(B1):2051. [4] Bak P, Tang C (1989) Earthquakes as a Self-Organized [11] MorenoY,CorreigAM,GomezJB,PachecoAF(2001) Critical Phenomenon. J Geophys Res 94:15635-15637. A model for complex aftershock sequences. J Geophys [5] Carlson J M, Langer J S. , Shaw B-E (1994) Dynamics Res B: Solid Earth 106:6609-6619. of earthquakefaults. Rev Mod Phys 66:657-670. [12] Helmstetter A., Shaw B. E. (2006) Relation between [6] BurridgeR,KnopoffL(1967)Modelandtheoreticalseis- stress heterogeneity and aftershock rate in the rate-and- micity.Bull Seismol Soc Am 57:341-362. state model. J Geophys Res 111:B07304. [7] OmoriF(1894)Ontheaftershocksofearthquakes.JColl [13] DieterichJH(1979) Modelingofrockfriction.1.Exper- 8 imental Results and Constitutive Equations. J Geophys B 58:6353-6366. Res 84:2161-2168. [20] Le Doussal P, Middleton A, Wiese K J (2008) Statis- [14] Ruina A (1983) Slip Instability and State Variable Fric- tics of static avalanches in a random pinning landscape. tion Laws. J GeophysRes 88:10359-10370. arXiv:0803.1142v1. [15] Ohmura A, Kawamura H (2007) Rate- and state- [21] Jagla E A (2007) Strain localization driven by struc- dependentfrictionlawandstatisticalpropertiesofearth- tural relaxation in sheared amorphous solids. Phys Rev quakes.**Europhys Lett 77:690011-5. E 76:0461191-7. [16] Fisher D S (1998) Collective transport in random me- [22] MaroneC(1998)Theeffectofloadingrateonstaticfric- dia: from superconductors to earthquakes. Physics Re- tion and the rate of fault healing during the earthquake ports 301:113-150. cycle. Nature 391:69-72. [17] LacombeF,ZapperiS,HermannHJ(2001)Forcefluctu- [23] B. N. J. Persson, Sliding Friction, Physical Principles ation in a drivenelastic chain.Phys Rev B 63:1041041-7. and Applications, (Springer, Berlin, 2000). [18] RossoA,LeDoussalP,WieseKJ(2007) Numericalcal- [24] By this, we mean all earthquakes in the region between culation of the functional renormalization group fixed- 115 and 119 W, and between 32 and 35 N, of M > 2 point functions at the depinning transition. Phys Rev B between January 1st (1980), and May 20th (2008), as 75:220201(R). reported at http://quake.geo.berkeley.edu/anss. [19] Zapperi S, Cizeau P, Durin G, Stanley H E (1998) Dy- [25] Marone C (1998) Laboratory-Derived Friction Laws and namics of a ferromagnetic domain wall: Avalanches, de- their Application to Seismic Faulting. Annu Rev Earth pinning transition, and the Barkhausen effect. Phys Rev Planet Sci 26:643-696.