Table Of ContentThe 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: jagla@cab.cnea.gov.ar 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.