Slip–weakening law

1. NMESD 2003 A FINITE DIFFERENCE ALGORITHM TO MODEL A FULLY 3–D DYNAMIC RUPTURE GOVERNED BY DIFFERENT CONSTITUTIVE LAWS Andrea Bizzarri, Massimo Cocco bizzarri@bo.ingv.it Istituto Nazionale di Geofisica e Vulcanologia, Italy Coupling of two modes of propagation Abstract In the last decades large effort has been expended to realize numerical algorithms for solving the fundamental elastodynamic equation to model earthquake ruptures. Two class of methods have been extensively implemented and used: the Boundary Integral Equation (BIE) and the Finite Difference (FD) approaches. In Bizzarri et al.(2001) we have compared 2–D pure in–plane solutions to discuss the main differences existing between these two different numerical strategies. In this work we present a new FD numerical method to solve the fully dynamic spontaneous problem for a truly 3–D rupture on planar faults. We implement the Traction–at–Split–Nodes of Andrews (1999) fault boundary condition for a system of faults, either vertical or oblique. For each fault we can assume different constitutive laws: we can use a slip–weakening law to prescribe the traction evolution within the breakdown zone or a rate– and state–dependent friction law, which involve the choice of a governing relation for the state variable. We have implemented a 3–D version of the Rosembrock Stiff Integration, generalizing that previously used in our 2–D model, and we have included the free surface effect. In our numerical procedure the initial shear stress is not necessarily imposed in only one spatial direction and the two components of slip, slip velocity and traction can vary during the dynamic crack propagation. Such components are coupled together in order to satisfy, in norm, the adopted governing law, and therefore they allow for dynamically controlled rake variations. We can also model faults with spatially heterogeneous distributions of constitutive parameters in order to simulate crack arrest and the healing of slip. At each time level the two components of slip, slip velocity and traction are calculated independently, with the requirements that their modula satisfy the assumed governing law and that the traction is collinear with the slip velocity. We did numerical simulations with initial rake of 45°, either with the slip–weakening and the Dieterich–Ruina law. We show that during the crack propagation the rake (i. e. the slip velocity azimuth) can vary, both in space (especially in the cohesive zone, after the tip bifurcation) and in time. As observed by Andrews (1994) for a mixed–mode crack, the variation of rake is lower as the initial value of stress increases, being the strength parameter S constant. All the simulations shown in the Figures refer to a single vertical fault. Frictional parameters Initial rake = 45 degrees Slip–weakening law with S = 0.8 Test #26 t0 = 1.0 tu = 1.8 tf = 0.0 d0 = 1.3 Test #30 t0 = 3.0 tu = 3.8 tf = 2.0 d0 = 1.3 Test #34 t0 = 4.0 tu = 4.8 tf = 3.0 d0 = 1.3 Test #35 t0 = 10.0 tu = 10.8 tf = 9.0 d0 = 1.3 Dieterich–Ruina law Test #11 a = 0.0085 b = 0.016 L = 10 mm The numerical method Our method allows for the calculation and interactions between a system of six faults, two vertical and four dipping faults. A grid of nodes is introduced by using quadrilateral isoparametric elements and the solutions of the fundamental elastodynamic equation is obtained by stepping explicitely the solutions at the previous time level through time: from displacement components at nodes ui we find strain tensor components at integration points, then stress tensor components and then net force components fi(n) at nodes:v(n)v(n– 1)a(n)Dt and u(n)u(n– 1)v(n)Dt, where a is the node acceleration derived from the second law of mechanics: a(n)f(n)/m, m being the mass matrix element. Such a formulation is equivalent to the local stiffness matrix. All components of slip velocity, displacement and force are defined at the same node point and the number of finite difference expressions that are calculated is 8 times larger than in a simple staggered grid scheme, in which velocities and stresses are defined at different points in space. In the latter way is easy to increase the order of accuracy to fourth order in space, while keeping it second order in time. The mass matrix is diagonal and the small displacement approximation is used: in such a way we can refer to the numerical algorithm either as a Finite Element (FE) scheme and Finite Difference (FD) scheme. For fault node point solutions are obtained by introducing the Traction–at–Split–Nodes (TSN) fault boundary condition which assumes an explicit displacement discontinuity between the fault surface. The fault strengthSfault is specified assuming a constitutive law (Sfaultm(Du, Dv, Y, …)sneff). For the linear slip–weakening law (Ida, 1972) we first find a trial value of traction on the fault plane that gives no differential acceleration (aaab, or Da 0):tt(n)t0 (mafa(n) – mbfb(n))/ (A ma + A mb), where A is the slip node area and maand mb are the split node masses, above and under the fault plane, respectively. If the fracture criterion is not satisfied, i. e. if Dv(n–½)vavb 0 tt(n) < Sfault then the fault slip Du is not changed and the traction t is not updated and the fault is not slipping in this case. On the contrary, if the fracture criterion is satisfied, the fault starts to slip or continues to slip. In this case we update the actual total dynamic traction by imposing the collinearity between traction and slip velocity vectors:t(n)c v*(n) where v*(n)Dv(n–½) + ½ ADt (1/ma + 1/mb) tt(n). From this equation we determine the director cosines: gi(n)vi*(n)/v*(n). Therefore the two actual components of total traction are expressed as: ti (n)gi(n)Sfault. For the whole class of rate– and state–dependent friction laws (Dieterich, 1986; Ruina, 1983; Roy and Marone, 1996) we solve for Dvi and Y the ODE system: a1 = daccdtau * (L1 – m(Dv, Y)(Dv1/Dv) sneff) a2 = daccdtau * (L2 – m(Dv, Y)(Dv2/Dv) sneff ) (d/dt)Y = evolution equation where daccdtau is the differential acceleration for unit stress drop and Li are traction that will give no change of differential velocity. We generalize to 3–D the Rosembrok stiff integration (Press et al., 1992), that perform better than Runge–Kutta algorithm. From Dvi(n) we calculate the director cosines gi(n)Dvi(n)/Dv(n) and finally the actual value of shear traction: ti (n)gi(n)Sfault. Vertical lines represent the time duration of the breakdown process: the first line indicates the time step at which the traction is at the upper yield value and the second one marks the time step at which it is at the kinetic level. Location #1 is along a line passing through the hypocenter and in the x direction; Location #2 is along a line at 45 degrees and Location #3 is along a line in the z direction. Effects of heterogeneities Dieterich–Ruina law Slip–weakening law Our methods allows for the modeling of heterogeneous configurations, in terms of initial stress distribution (pre–stress) and in terms of constitutive parameters. We present here two numerical experiments, one with the slip–weakening law and a barrier (a region with a high value of strength parameter S) and one with Dieterich–Ruina law and a velocity strengthening region, surrounding a velocity weakening one. In such cases the coupling of the two modes of propagation is quite different with respect to the homogeneous configurations and the rake variation is quite complex, especially in the case of crack arrest. a = 0.012 b = 0.016 S = 19 a = 0.014 b = 0.012 Reference cases S = 0.8 We compare solutions obtained with our 2–D pure in–plane FD code and with this new fully 3–D FD algorithm, adopting either a linear slip–weakening law and a Dieterich–Ruina law. In this simulation, for the 3–D model the initial stress is non null along the in–plane direction only. Nucleation strategies are the same for two methods (Bizzarri et al., 2001). Our numerical experiments show that the bifurcation phenomenon is quite different and the superpositions of the solutions along the time generalize the results of Tinti et al. (2003). Conclusions We have proposed a new Finite Difference numerical method to model the truly 3–D (not mixed–mode) fully dynamic and spontaneous crack propagation on a system of six mutually interacting faults, embedded in an elastic medium, with homogeneous or heterogeneous properties. The rheology on fault planes is described by different friction laws: the linear slip–weakening and the rate– and state – dependent friction laws. Our methodology allows for the rake rotation during the rupture propagation. We have shown that such spatial and temporal variations occour within the breakdown zone, where the physical processes of the traction drop and the seismic wave emission take place, especially when the crack tip bifurcates. We have also shown that the behavior of the solution is quite different adopting a linear slip–weakening constitutive equation or a Dieterich–Ruina law, generalizing the results of Bizzarri et al. ( 2001). We are able to simulate real–world events, accounting for the spatially heterogeneous initial stress, the free surface effect and the heterogeneities of the frictional parameters on the fault planes. • Model parameters • Slip–weakening law (adimensional set) • = 1. vS = 1. vP = 1.732 Dx = Dx = Dx = 0.2 t0 = 1.0 tu = 1.8 tf = 1.0 d0 = 1.3 • Dieterich–Ruina law • = 2700 Kg/m3vS = 3000 m/svP = 5196 m/sDx = Dx = Dx = 0.1 m • a = 0.008 b = 0.016 L = 10 mm sneff = 100 MPa • vinit = 10 mm/s References Andrews D. J. (1994): Dynamic growth of mixed–mode shear cracks, Bull. Seism. Soc. Am., 84, No. 4, pp. 1184–1198. Andrews D. J. (1999): Test of two methods for faulting in finitedifference calculations, Bull. Seism. Soc. Am., 89, No. 4, pp. 931937. Bizzarri A., Cocco M., Andrews D. J., Boschi E. (2001): Solving the dynamic rupture problem with different numerical approaches and constitutive laws, Geophys. J. Int., 144, pp. 656–678. Dieterich J. H. (1986): A model for the nucleation of earthquake slip, Earthquake Source Mechanics, Geophysical Monograph, 37, S. Das, J. Boatwright and C. H. Scholz eds., AGU, Washington, DC, pp. 37–47. Ida Y. (1972): Cohesive force across the tip of a longitudinalshear crack and Griffith’ s specific surface energy, J. Geophys. Res., 77, No. 20, pp. 3796–3805. Press W. H., Teukolsky S. A., Vetterling W. T., Flannery B. P. (1992): Numerical Recipes in FORTRAN. The art of scientific computing, Second Edition, Cambridge University Press. Roy M., Marone C. (1996): Earthquake nucleation on model faults with rate– and state–dependent friction: effects of inertia, J. Geophys. Res., 101, No. B6, pp. 13919–13932. Ruina A. L. (1983): Slip instability and state variable friction laws, J. Geophys. Res., 88, No. B12, pp. 10359–10370.