VIEWS: 3 PAGES: 30 POSTED ON: 2/24/2012
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 2. Aftershock decay rates can be the mainshock rupture zone 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 4. The largest aftershock magnitude is typically about 1-1.5 units below 5. Aftershocks behavior is that of the mainshock (Båth law). NOT universal! 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. Non-linear Continuum Damage Rheology (1) Mechanical aspect: sensitivity of elastic moduli to distributed cracks and sense of loading. peak stress yielding Stress a=0 Strain 0 < a < ac s s Tension Tension Compression e 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 I1 U = I 1 I 2 - I 1 I 2 = 2 I2 Where and are Lame constants; I1= ekk is an additional elastic modulus I2= eijeij U I2 I1 sij = = - I1ij 2 - eij eij I1 I2 Non-linear Continuum Damage Rheology (2) Kinetic aspect associated with damage evolution peak stress yielding Stress Strain 0 < a < ac 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 = F TS = s ij e ij - i J i dU d 1 Energy balance dt dt Entropy balance dS Ji = - i G dt T Gibbs equation 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 T a dt da da/dt > 0 = Cd I 2 - 0 dt > 0 I1 Strain weakening = invariant Shear I2 ratio (degradation) Stress = tan () sn = 0 I1=ekk I2= eijeij da/dt < 0 da a < 0 = C1 exp( ) I 2 - 0 dt C2 healing (strengthening) Normal Stress sn Rate- and state-dependent friction experiments constrain parameters c1 and c2. D 0.1 1 10 V / V0 D = D/s For details, see Lyakhovsky 0 0.1 0.2 0.3 0.4 0.5 et al. (GJI, 2005) LOGe (s/s0) 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 500 G3 data Simulation 400 Y Stress, MPa 300 1 = , a 0 X 200 Cva 100 1/Cv) = 5·1010 Pa, Cd = 3 s-1 0 0 2 4 6 8 10 12 14 Strain ( x103 ) Berea sandstone under 50 MPa confining pressure 250 200 Differential stress (MPa) 150 100 50 Accumulated irreversible strain 0 -1 -0.5 0 0.5 1 1.5 0 = 1.4 1010 Pa, Strain % Cv = 10-10 Pa-1, 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 = 20(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-to- failure 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, 2de/dt = –e (4) de da Using this condition in (2) and (3) gives = -C v 0 1 - a e (5) dt dt 1 2 and integrating (5) we get e = A exp R1 - a (6) 2 1 2 where R = d/M = 0Cv and A = e s exp - R1 - 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 dt = C d e s2 exp R1 - a - R1 - a s - e 0 2 2 2 (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 f dN dt = Cd e s2 exp R1 - a s - fN - R1 - a s - e 0 2 2 2 (9) If fN is small (generally true), so that (fN)2 can be neglected f dN dt = Cd e s2 exp- 2fNR1 - a s - e 0 2 (10) If also the initial strain induced by the mainshock is large enough so that e 02 << e s2 exp - 2fNR 1 - a s (11) the solution is (the Omori-Utsu law) dN Cd e s2 = dt 2fR1 - a s Cd e s2t f (12) C d e s2 For t = 0 N0 = f dN N0 N0 1 so = = (13) t 1 2fR1 - a N t 1 2fR1 - a N dt 2fR1 - a s N 0 s 0 s 0 The parameters of the Omori-Utsu law are k= 1 dN/dt = K(c + t)-p 2fR1 - a s 1 k c= = 2fR1 - a s N 0 N 0 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 dt = Cd e s2 exp R1 - a s - fN - R1 - a s - e 0 2 2 2 (9) Events rate vs. time for several values of R = d/M with e0=0, as=0) Timescale of fracturing Material property R = Timescale of stress relaxation 180 Small R: 160 •expect long R = 0.1 active Number of aftershocks per day 140 aftershock 120 sequences 100 R=1 80 Modified Omori 60 law with p = 1 Large R: 40 •expect short R = 10 diffuse 20 sequences 0 0 10 20 30 40 50 60 70 80 90 100 Time (day) 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) Timescale of fracturing Material property R = Timescale of stress relaxation 150 Omori p=1 125 Number of aftershocks per day 100 R = 0.1 Omori 75 p=1 R = 0.3 50 Omori p = 1.2 25 R = 10 R=1 0 0 20 40 60 80 100 Time (day) Imposed damage 3-D numerical simulations (major fault zone) 1-7 km 35 km Newtonian Sedimentary cover viscosity Damage visco-elastic rheology Crystalline plus power law viscosity Crust (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 400 500 0 Initial regional stress for temperature gradients 20 oC/km – heavy line 10 s = s n 30 oC/km – dash line 40 oC/km – dotted line Strain rate = 10-15 1/s 20 Brittle-ductile e = A0 e -Q / RT n transition at 300 oC 30 Moho 40 50 Simulations with fixed r = 300 s Effects of R (sediment thickness = 1 km, gradient 20 oC/km ) 200 200 100 R=0.1 R=1 Number of events Number of events R=2 Number of events 150 80 150 60 100 100 40 50 50 20 0 0 0 0 20 40 60 80 100 4060 2080 0 100 0 20 40 60 80 100 Time (days) Time (days) Time (days) 50 50 Increasing R values: 45 45 40 R=3 40 Number of events R=10 Number of events 35 35 30 25 30 25 •diffuse sequences 20 20 15 15 •shorter duration 10 10 5 5 0 0 •smaller # of events 0 20 40 60 80 100 0 20 40 60 80 100 Time (days) Time (days) 4 R=0.1 3.5 R=1 R=2 3 R=3 2.5 Log(Number) 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 20 oC/km ) 200 250 100 S=1 km S=4 km S=7 km Number of events Number of events Number of events 200 80 150 150 60 100 100 40 50 50 20 0 0 0 0 20 40 60 80 100 0 20 40 60 80 100 0 20 40 60 80 100 Time (days) Time (days) 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 -10 Depth (km) -15 -20 -25 -30 R = 0.1 -35 T = 20 C/km -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) “Regional” depth (1283 events) d95 d5% JV HypoDD Johnson Valley Fault 11944 events 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 s = 2 e seis can be estimated as The rate of gradual inelastic strain can de i dt = -aCvs / 2 be estimated as The inelastic strain accumulation (aseismic creep) is e i = Cvs / 2 R Slip ratio Seismic slip e seis 1 = = 0.1 90 % Total slip e total 1 R 1 50 % 2 33 % 10 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, 1155- 1165, 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.