A generalized law for aftershock behavior in a damage rheology model

1. A generalized law for aftershock behavior in a damage rheology model • 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 Yehuda Ben-Zion1 and Vladimir Lyakhovsky2 1. University of Southern California 2. Geological Survey of Israel

2. Main observed features of aftershock sequences: 1. Aftershocks occur around the mainshock rupture zone 2. Aftershock decay rates can be described bythe 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 that of the mainshock (Båth law). 5.Aftershocks behavior is NOT universal!

3. 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.

4. peak stress yielding Stress 0 < a < ac a = 0 Strain s s Tension Tension e e Compression Compression Non-linear Continuum Damage Rheology(1)Mechanical aspect: sensitivity of elastic moduli to distributed cracks and sense of loading.

5. This is accounted for by generalizing the strain energy function of a deforming solid The elastic energy U is written as: I1= kk I2= ijij Where  and  are Lame constants;  is an additional elastic modulus

6. peak stress yielding Stress a = ac Strain s s Tension Tension e e Compression Compression Non-linear Continuum Damage Rheology(2)Kinetic aspect associated with damage evolution 0 < a < ac

7. This is accounted for by making the elastic modulifunctions of a damage state variable(x, y, z, t), representing crack density in a unit volume, and deriving an evolution equation fora.

8. 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 The internal entropy production rate per unit mass, G, is:

9. Strain invariant ratio I1=kk I2= ijij

10. Rate- and state-dependent friction experiments constrain parameters c1 and c2. For details, see Lyakhovsky et al. (GJI, 2005)

11. 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) The Maxwell relaxation time is as small as tens of years Viscosity = 8 x 1019 Pa s

12. Non-linear Continuum Damage Rheology(3)damage-related viscosity Stress-strain and AE locations for G3 (Lockner et al., 1992) Z Y X (1/Cv) = 5·1010 Pa, Cd = 3 s-1

13. Berea sandstone under 50 MPa confining pressure 250 200 150 Differential stress (MPa) 100 Accumulated irreversible strain 50 0 -1 -0.5 0 0.5 1 1.5 m0 = 1.4 1010 Pa, Cv = 10-10 Pa-1, R = 1.4 Strain % Data from Lockner lab. USGS Model from Hamiel et al., 2004

14. What aboutaftershocks?

15. Aftershocks: 1D analytical results for uniform deformation For 1D deformation, the equation for positive damage evolution is da/dt = Cd (e2-e02), (1) whereeis the current strain and e0 separates degradation from healing. The stress-strain relation in this case is s = 2m0(1 – a)e, (2) where m0(1–a) is the effective elastic modulus of a 1D damaged material with m0 being the initial modulus of the undamaged solid. (Ben-Zion and Lyakhovsky [2002] showed analytically that these equations lead underconstant stressloading 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= (Cvda/dt)s(3)

16. For aftershocks, we consider material relaxation following a strain step. This corresponds to a situation with a boundary conditions ofconstant total strain. In this case the rate of elastic strain relaxation is equal to the viscous strain rate, 2de/dt = –e(4) Using this condition in (2) and (3) gives (5) and integrating (5) we get (6) whereR =td/tM=m0Cvand is integration constant with a =as and e =es for t = 0. Using these results in (1) yieldsexponentialdamage evolution (7)

17. Scaling the results to number of events N Assuming that a is scaled linearly with the number of aftershocks N (8) we get (9) If fN is small (generally true), so that (fN)2 can be neglected (10) If also the initial strain induced by the mainshock is large enough so that (11) the solution is (the Omori-Utsu law) (12)

18. For t = 0 so (13) The parameters of the Omori-Utsu law are and p = 1 We now return to the general exponential equation(9) and examine analytical results first withe0=0, as=0 and then with finite small values. (9) dN/dt = K(c + t)-p

19. Events rate vs. time for several values ofR =td/tM (with e0=0, as=0) 180 160 • Small R: • expect long active aftershock sequences 140 120 Number of aftershocks per day 100 80 Modified Omori law with p = 1 60 40 • Large R: • expect short diffuse sequences 20 0 0 10 20 30 40 50 60 70 80 90 100 Time (day) Timescale of fracturing Material propertyR = Timescale of stress relaxation R = 0.1 R = 1 R = 10 Changing thepower-lawparameters,we canfit theother lines !!!

20. 150 125 100 Number of aftershocks per day 75 50 25 0 0 20 40 60 80 100 Time (day) Events rate vs. time for several values ofR =td/tM (with finite e0,as) Timescale of fracturing Material propertyR = Timescale of stress relaxation Omori p = 1 R = 0.1 Omori p = 1 R = 0.3 Omori p = 1.2 R = 1 R = 10

21. Imposed damage (major fault zone) 3-D numerical simulations 1-7 km Free surface Basement 35 km Newtonian viscosity Sedimentary cover Moho Damage visco-elastic rheology plus power law viscosity (based on diabase lab data) Crystalline Crust 50 km Upper mantle 100 km Damage visco-elastic rheology plus power-law viscosity (based on Olivine lab data) y 100 km x z In each layer the strain is the sum of damage-elastic, damage-related inelastic, and ductile components: Initial stress = regional stress + imposed mainshock slip on a fault extending over 50 km ≤ y ≤ 150 km, 0 ≤ z ≤ 15 kmwith fixed boundaries

22. 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 Brittle-ductile transition at 300 oC

23. Simulations with fixedtr = 300 s Effects of R (sediment thickness = 1 km, gradient 20 oC/km) 200 200 100 R=0.1 R=1 R=2 80 150 150 60 100 Number of events 100 Number of events Number of events 40 50 50 20 0 0 0 40 60 80 100 0 20 40 60 80 100 0 20 40 0 20 60 80 100 Time (days) Time (days) Time (days) 50 50 45 45 • Increasing R values: • diffuse sequences • shorter duration • smaller # of events R=3 R=10 40 40 35 35 30 30 25 25 Number of events Number of events 20 20 15 15 10 10 5 5 0 0 0 20 40 60 80 100 0 20 40 60 80 100 Time (days) Time (days)

24. 4 R=0.1 R=1 3.5 R=2 3 R=3 2.5 R=10 2 Log(Number) 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

25. Effect of Sediment thickness(R = 1, gradient 20 oC/km) 200 250 100 S=1 km S=4 km S=7 km 200 80 150 150 60 Number of events Number of events Number of events 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)

26. 0 -5 -10 -15 -20 Depth (km) -25 -30 R = 0.1 T = 20 C/km -35 -40 0 30 60 90 120 150 180 210 240 270 300 Time (days) Effect of thermal gradient and R (sediment layer 1 km) Increasing thermal gradient and/or R: thinner seismogenic zone The maximum event depth decreases with time from the mainshock

27. HypoDD Hauksson (2000) JV Depth of seismic-aseismic transition increases following Landers EQ and then shallows by ≤ 3 km over the course of 4 yrs. Observed Depth Evolution of Landers aftershocks (Rolandone et al., 2004) d95 “Regional” depth (1283 events) d5% Johnson Valley Fault 11944 events

28. 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 The inelastic strainaccumulation(aseismic creep) is Seismic slip Total slip

29. 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 ratioR of the timescale for brittlefracture 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 powerlaw 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.

30. Thankyou 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.