Transcript Spatio-Temporal Evolution of Earthquakes and Faults

A generalized law for aftershock
behavior in a damage rheology model
Yehuda Ben-Zion1 and Vladimir Lyakhovsky2
1. University of Southern California
2. Geological Survey of Israel
•
•
•
•
•
Outline
Brief background on aftershocks
Brief background on the employed damage rheology
1-D Analytical results on aftershocks
3-D Numerical results on aftershocks
Discussion and Conclusions
Main observed features of aftershock sequences:
1. Aftershocks occur around
the mainshock rupture zone
2. Aftershock decay rates can be
described by the Omori-Utsu law:
DN/Dt = K(c + t)-p
However, aftershock decay rates can
also be fitted with exponential and
other functions (e.g., Kisslinger, 1996).
3. The frequency-size statistics of
aftershocks follow the GR relation:
logN(M) = a - bM
5. Aftershocks behavior is
NOT universal!
4. The largest aftershock magnitude
is typically about 1-1.5 units below
that of the mainshock (Båth law).
Existing aftershock models:
•Migration of pore fluids (e.g., Nur and Booker, 1972)
•Stress corrosion (e.g., Yamashita and Knopoff, 1987)
•Criticality (e.g., Bak et al., 1987; Amit et al., 2005)
•Rate- and state-dependent friction (Dieterich, 1994)
•Fault patches governed by dislocation creep (Zöller et al., 2005).
Is the problem solved?
The above models focus primarily on rates.
None explains properties (1)-(5), including the observed spatio-temporal
variability, in terms of basic geological and physical properties.
This is done here with a damage rheology framework and realistic model
of the lithosphere.
Stress
Non-linear Continuum Damage Rheology (1) Mechanical aspect:
sensitivity of elastic moduli to distributed cracks and sense of
loading.
peak
stress
yielding
a=0
Strain
s
0 < a < ac
Tension
Compression
e
s
Tension
Compression
e
This is accounted for by generalizing the strain
energy function of a deforming solid
The elastic energy U is written as:
1 2

U =  I 1  I 2 - I 1 I 2 

 2
Where  and  are Lame constants;
 is an additional elastic modulus
I1= ekk
I2= eijeij
I1
=
I2


I2 
U 
I
1
 I1ij   2 - 
eij
sij = 
= -


e ij 
I1 
I
2 

Stress
Non-linear Continuum Damage Rheology (2) Kinetic aspect associated
with damage evolution
peak
stress
yielding
0 < a < ac
Strain
a = ac
s
s
Tension
Tension
Compression
e
Compression
e
This is accounted for by making the
elastic moduli functions of a damage
state variable a(x, y, z, t), representing
crack density in a unit volume, and
deriving an evolution equation for a.
Thermodynamics
Free energy of a solid, F, is
F = F(T, eij, a)
T – temperature, eij – elastic strain tensor,
a – scalar damage parameter
Energy balance
Entropy balance
Gibbs equation
dU d
1
= F  TS = s ij e ij -  i J i
dt dt

dS
 Ji 
= - i    G
dt
T
F
F
dF = -SdT 
de ij 
da
e ij
a
The internal entropy production rate per unit mass, G, is:
Ji
1
1 F da
G = - 2  i T  s ij e ij 0
T
T a dt
T
da
= Cd  I 2  - 0 
dt
da/dt > 0
 > 0
Shear
Stress

I1
=
I2
weakening
(degradation)
 = tan () sn
 = 0
da/dt < 0
 < 0
healing
(strengthening)
Normal Stress
sn
Strain
invariant
ratio
I1=ekk
I2= eijeij
da
a
= C1  exp( )  I 2  - 0 
dt
C2
D
D = D/s
Rate- and state-dependent
friction experiments constrain
parameters c1 and c2.
For details, see Lyakhovsky
et al. (GJI, 2005)
ss
Non-linear Continuum Damage Rheology (3) damage-related viscosity
10 years creep experiment
on Granite beam at room
temperature
Ito & Kumagai, 1994
For typical values of
shear moduli of granite
(2-3 * 1010 Pa)
Viscosity = 8 x 1019 Pa s
The Maxwell relaxation time
is as small as
tens of years
Non-linear Continuum Damage Rheology (3) damage-related viscosity
Stress-strain and AE locations
for G3 (Lockner et al., 1992)
600
Z
G3 data
Simulation
500
X
Stress, MPa
Y
400
300
1
=
, a  0
Cva
200
1/Cv) = 5·1010 Pa, Cd = 3 s-1
100
0
0
2
4
6
8
Strain ( x103 )
10
12
14
Berea sandstone under 50 MPa confining pressure
Differential stress (MPa)
250
200
150
100
Accumulated
irreversible
strain
50
0
-1
0 = 1.4 1010 Pa,
Cv = 10-10 Pa-1,
-0.5
0
0.5
1
1.5
Strain %
Data from Lockner lab. USGS
Model from Hamiel et al., 2004
R = 1.4
What about aftershocks?
Aftershocks: 1D analytical results for uniform deformation
For 1D deformation, the equation for positive damage evolution is
da/dt = Cd (e2-e02),
(1)
where e is the current strain and e0 separates degradation from healing.
The stress-strain relation in this case is
s = 20(1 – a)e,
(2)
where 0(1–a) is the effective elastic modulus of a 1D damaged material
with 0 being the initial modulus of the undamaged solid.
(Ben-Zion and Lyakhovsky [2002] showed analytically that these
equations lead under constant stress loading to a power law time-tofailure relation with exponent 1/3 for a system-size brittle event).
For positive rate of damage evolution (e > e0), we assume inelastic strain
before macroscopic failure in the form
e = (Cv da/dt) s
(3)
For aftershocks, we consider material relaxation following a strain step.
This corresponds to a situation with a boundary conditions of
constant total strain.
In this case the rate of elastic strain relaxation is equal to the viscous
strain rate,
(4)
2de/dt = –e
Using this condition in (2) and (3) gives
and integrating (5) we get
where R = d/M = 0Cv and
de
da
= -C v  0 1 - a   e
dt
dt
1
2
e = A  exp R1 - a  
2

(5)
(6)
 1
2
A = e s  exp- R1 - a s   is integration
 2

constant with a = as and e = es for t = 0.
Using these results in (1) yields exponential damage evolution


 
da
2
2
= C d  e s2 exp R1 - a  - R1 - a s  - e 02
dt
(7)
Scaling the results to number of events N
Assuming that a is scaled linearly with the number of aftershocks N
a = a s  fN
(8)
we get


 
dN
2
2
f
= Cd  e s2 exp R1 - a s - fN  - R1 - a s  - e 02
dt
(9)
If fN is small (generally true), so that (fN)2 can be neglected
f

dN
= Cd  e s2 exp - 2fNR 1 - a s  - e 02
dt

(10)
If also the initial strain induced by the mainshock is large enough so that
e 02 << e s2 exp- 2fNR1 - a s 
the solution is (the Omori-Utsu law)
Cd e s2
dN
=
dt 2fR1 - a s Cd e s2t  f
(11)
(12)
N 0 =
For t = 0
so
Cd e s2
f
N 0
N 0
dN
1
=
=



dt 2fR1 - a s N 0t  1 2fR1 - a s N 0 t  1 2fR1 - a s N 0
(13)
The parameters of the Omori-Utsu law are
dN/dt = K(c + t)-p
1
2fR1 - a s 
1
k
c=
=
2fR1 - a s N 0 N 0
k=
and p = 1
We now return to the general exponential equation (9) and examine
analytical results first with e0=0, as=0 and then with finite small values.
f


 
dN
2
2
= Cd  e s2 exp R1 - a s - fN  - R1 - a s  - e 02
dt
(9)
Events rate vs. time for several values of R = d/M with e0=0, as=0)
Material property R =
Timescale of fracturing
Timescale of stress relaxation
180
Small R:
•expect long
active
aftershock
sequences
Number of aftershocks per day
160
Large R:
•expect short
diffuse
sequences
R = 0.1
140
120
100
R=1
80
Modified Omori
law with p = 1
60
40
R = 10
20
0
0
10
20
30
40
50
Time (day)
60
70
80
90
100
Changing the power-law parameters, we can fit the other lines !!!
Events rate vs. time for several values of R = d/M with finite e0,as)
Material property R =
Timescale of fracturing
Timescale of stress relaxation
150
Number of aftershocks per day
Omori
p=1
125
R = 0.1
100
Omori
p=1
75
50
R = 0.3
Omori
p = 1.2
25
R=1
R = 10
0
0
20
40
60
Time (day)
80
100
Imposed damage
(major fault zone)
3-D numerical simulations
1-7 km
35 km
Newtonian
viscosity
Sedimentary cover
Crystalline
Crust
Damage visco-elastic rheology
plus power law viscosity
(based on diabase lab data)
50 km
Upper mantle
Damage visco-elastic rheology
plus power-law viscosity
(based on Olivine lab data)
100 km
y
x
z
In each layer the strain is the sum of damage-elastic, damage-related
inelastic, and ductile components: e t = e e  e i  e d
ij
ij
ij
ij
Initial stress = regional stress + imposed mainshock slip on a fault extending
over 50 km ≤ y ≤ 150 km, 0 ≤ z ≤ 15 km with fixed boundaries
Differential Stress (MPa)
0
100
200
300
0
 s = s n
10
400
500
Initial regional stress for
temperature gradients
20 oC/km – heavy line
30 oC/km – dash line
40 oC/km – dotted line
Strain rate = 10-15 1/s
20
e = A0e

-Q / RT n
30
Moho
40
50
Brittle-ductile
transition at 300 oC
Simulations with fixed r = 300 s
Effects of R (sediment thickness = 1 km, gradient
200
100
50
0
0
20
40
60
80
100
R=1
150
100
50
0
4060
0
2080
100
50
45
40
35
30
25
20
15
10
5
0
R=3
20
100
0
20
40
Time (days)
40
60
Time (days)
80
Number of events
Number of events
Time (days)
0
Number of events
R=0.1
150
Number of events
Number of events
200
100
20 oC/km )
50
45
40
35
30
25
20
15
10
5
0
R=2
80
60
40
20
0
60
80
Time (days)
Increasing R values:
R=10
•diffuse sequences
•shorter duration
0
20
40
60
Time (days)
80
100
•smaller # of events
100
4
R=0.1
R=1
3.5
R=2
3
Log(Number)
R=3
2.5
R=10
2
1.5
1
0.5
0
3
4
5
6
Magnitude
Small R values (R < 1): Power law frequency-size statistics
Large R values (R > 3): Narrow range of event sizes
Effect of Sediment thickness (R = 1, gradient
250
100
S=4 km
200
Number of events
S=1 km
150
Number of events
Number of events
200
150
100
100
50
0
50
0
0
20
40
60
Time (days)
80
100
20 oC/km )
S=7 km
80
60
40
20
0
0
20
40
60
Time (days)
80
100
0
20
40
60
80
100
Time (days)
Increasing thickness of weak sediments: diffuse sequences, shorter
duration, smaller number of events (similar to increasing R values)
Effect of thermal gradient and R (sediment layer 1
km )
0
-5
Depth (km)
-10
-15
-20
-25
-30
R = 0.1
T = 20 C/km
-35
-40
0
30
60
90
120
150
180 210
240 270
300
Time (days)
Increasing thermal gradient and/or R: thinner seismogenic zone
The maximum event depth decreases with time from the mainshock
Observed Depth Evolution of Landers aftershocks
(Rolandone et al., 2004)
d95
d5%
“Regional” depth
(1283 events)
JV
Johnson
Valley
11944
eventsFault
HypoDD
Hauksson (2000)
Depth of seismic-aseismic transition increases following Landers EQ
and then shallows by ≤ 3 km over the course of 4 yrs.
The parameter R controls the partition of energy between seismic
and aseismic components (degree of seismic coupling across a fault)
The brittle (seismic) component of deformation
can be estimated as
The rate of gradual inelastic strain can
be estimated as
de i dt = -aCvs / 2
The inelastic strain accumulation (aseismic creep) is
Seismic slip
Total slip
e seis
1
=
=
e total 1  R
s = 2  e seis
e i = Cvs / 2
R
0.1
1
2
10
Slip ratio
90 %
50 %
33 %
10%
Main Conclusions
•Aftershocks decay rate may be governed by exponential rather than power
law as is commonly believed (see also Dieterich, 1994; Gross and Kisslinger,
1994, Narteau et al., 2002)
•The key factor controlling aftershocks behavior is the ratio R of the
timescale for brittle fracture evolution to viscous relaxation timescale.
•The material parameter R increases with increasing heat and fluids, and is
inversely proportional to the degree of seismic coupling.
•Situations with R ≤ 1, representing highly brittle cases, produce clear
aftershock sequences that can be fitted well by the Omori power law
relation with p ≈ 1, and have power law frequency size statistics.
•Situations with R >> 1, representing stable cases with low seismic coupling,
produce diffuse aftershock sequences & swarm-like behavior.
•Increasing thickness of weak sedimentary cover produce results that are
similar to those associated with increasing R.
Thank you
Key References (on damage and evolution of earthquakes & faults):
Lyakhovsky, V., Y. Ben-Zion and A. Agnon, Distributed Damage, Faulting, and Friction, J.
Geophys. Res., 102, 27635-27649, 1997.
Ben-Zion, Y., K. Dahmen, V. Lyakhovsky, D. Ertas and A. Agnon, Self-Driven Mode Switching
of Earthquake Activity on a Fault System, Earth Planet. Sci. Lett., 172/1-2, 11-21, 1999.
Lyakhovsky, V., Y. Ben-Zion and A. Agnon, Earthquake Cycle, Fault Zones, and Seismicity
Patterns in a Rheologically Layered Lithosphere, J. Geophys. Res., 106, 4103-4120, 2001.
Ben-Zion, Y. and V. Lyakhovsky, Accelerated Seismic Release and Related Aspects of
Seismicity Patterns on Earthquake Faults, Pure Appl. Geophys., 159, 2385 –2412, 2002.
Hamiel, Y., *Liu, Y., V. Lyakhovsky, Y. Ben-Zion and D. Lockner, A Visco-Elastic Damage
Model with Applications to Stable and Unstable fracturing, Geophys. J. Int., 159, 11551165, doi: 10.1111/j.1365-246X.2004.02452.x, 2004.
Ben-Zion, Y. and V. Lyakhovsky, Analysis of Aftershocks in a Lithospheric Model with
Seismogenic Zone Governed by Damage Rheology, Geophys. J. Int., in press, 2006.